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;  (4)  INTRODUCTION 

The  goal  of  the  “Breast  Cancer  Screening  Using  Photonic  Technology”  research  project  was  to  develop 
optical  imaging  techniques  that  make  use  of  noninvasive  near-infrared  (NIR)  light  for  obtaining  direct  two- 
dimensional  (2 -D),  and  tomographic  three-dimensional  (3-D)  images  of  cancerous  lesions  of  human  breast.  The 
imaging  method  involved  illuminating  the  specimen  with  ultrashort  NIR  pulses  of  laser  light  and  construction 
of  images  using  two  approaches.  The  time-sliced  direct  or  the  shadowgram  method  utilized  the  image  bearing 
component  of  the  forward-transmitted  light,  while  the  inverse  reconstruction  method  made  use  of  the 
measured  transmitted,  forward-scattered  or  backscattered  light  intensity  profiles,  known  or  estimated  optical 
properties  of  the  sample,  a  model  for  light  propagation  through  turbid  media  and  a  sophisticated  computer 
algorithm  to  construct  images  of  the  interior  structure  of  the  specimen. 

Significant  advances  were  made  not  only  in  developing  and  testing  both  of  these  approaches,  but  also  in 
demonstrating  the  diagnostic  potential  of  the  NIR  spectroscopic  imaging,  as  well  as,  deriving  analytical 
solutions  of  radiative  transport  equation  (RTE)  and  development  of  image  reconstruction  algorithms  for 
constructing  3-D  images.  In  order  to  assess  the  efficacy  of  the  photonics  approach,  the  photonics  experimental 
results  were  compared  and  correlated  with  that  from  histopathology  and  nuclear  magnetic  resonance  (NMR). 

The  results  are  intriguing  and  photonics  approach  highly  promising. 

This  final  report  covering  the  research  carried  out  during  the  entire  period  of  the  grant  is  organized  as  follows. 
An  overview  of  the  work  reported  in  the  three  Annual  Reports1'3  will  be  presented  for  completeness,  and  the 
details  will  be  referred  to  those  reports.  The  work  carried  out  since  the  last  report  will  be  presented  in  more 
detail.  References  to  appended  publications  will  be  made  for  detailed  descriptions  of  specific  technical  aspects 
of  the  research.  The  research  accomplishments  associated  with  tasks  outlined  in  the  approved  Statement  of 
Work  (SOW)  will  be  presented  in  the  “BODY”  Section.  “KEY  RESEARCH  ACCOMPLISHMENTS”  will 
present  a  bulleted  list  of  the  major  accomplishments  resulting  from  this  research.  A  list  of  publications, 
presentations,  employment  and  research  opportunities  and  other  deliverables  emanating  from  this  research  will 
appear  in  the  “REPORTABLE  OUTCOME”  Section.  The  importance  and  implications  of  the  completed 
research  will  be  summarized  in  the  “CONCLUSIONS”  Section.  Finally,  all  references  pertinent  to  the  report 
will  be  presented  in  the  “REFERENCES”  Section. 

(5)  BODY 

The  research  carried  out,  tasks  performed,  and  progress  made  over  the  duration  of  the  project  may  be  broadly 
grouped  as  follows: 

5.1.  Setting  up  and  Evaluation  of  Experimental  Arrangements; 

5.2.  Time-sliced  Direct  Imaging  of  ex  vivo  Normal  and  Cancerous  Human  Breast  Tissues; 

5.3.  Spectroscopic  Direct  Imaging  of  ex  vivo  Normal  and  Cancerous  Human  Breast  Tissues; 

5.4.  Development  of  Analytical  Solutions  of  Radiative  Transport  Equation; 

5.5.  Development  and  Demonstration  of  Inverse  Image  Reconstruction;  and 

5.6.  Correlation  with  Existing  Methods. 

We  will  present  a  brief  description  of  our  accomplishments  in  each  of  these  areas  and  refer  to  appended  articles 
and  reports  for  details. 

5.1 .  Setting  up  and  Evaluation  of  Experimental  Arrangements 

A  major  thrust  in  the  initial  phase  of  the  project,  as  detailed  in  the  approved  Statement  of  Work  (SOW),  was  to 
set  up,  and  evaluate  two  different  types  of  experimental  arrangements  for  time-sliced  imaging  [Technical 
Objectives  (TO)  1,2,  and  3;  Tasks  1-7].  Both  the  arrangements  used  a  femtosecond  Ti: sapphire  laser  and 
amplifier  system4  for  sample  illumination,  but  differed  in  the  signal  detection  scheme.  One  of  the  detection 
schemes  was  based  on  an  ultrafast  gated  intensified  camera  system  (UGICS),  which  provided  an  electronic  time 
gate  with  a  minimum  gate-width  of  80  ps  and  whose  position  could  be  varied  over  a  20-ns  interval  [T03,  Tasks 
5-7].  The  second  approach  [TOl,  Tasks  1-3]  used  an  optical  Kerr  gate  (OKG),5  along  with  a  Fourier  space 
gate,6  and  a  polarization  gate.7  The  setting  up,  testing,  and  evaluation  of  experimental  arrangements  [T02, 
Task  4]  using  both  the  schemes  were  detailed  in  the  First  Annual  Report.1  Our  comparison  [T03,  Task  7]  of 


4 


'  * 

•  the  efficacy  of  the  two  schemes  demonstrated  that  the  OKG-based  system  provided  higher  temporal  and  spatial 
resolution,  while  the  UGICS  system  was  more  user  friendly,  easier  to  implement,  and  provided  higher  signal 
level  since  the  gate  width  was  wider.  It  was  also  concluded  that  for  breast  imaging  the  resolution  provided  by 
the  UGICS  would  be  acceptable  given  the  competing  requirement  for  adequate  signal  level.  While  a  narrow 
time  gate  improves  the  resolution,  it  also  reduces  the  signal  level  that  is  crucial  for  forming  a  well-defined 
image.  Consequently,  we  mainly  used  the  UGICS-based  system  for  time-resolved  imaging  experiments. 

Another  important  experimental  approach  that  we  developed  and  used  in  the  project  was  the  spectroscopic 
imaging  arrangement,  also  detailed  in  the  First  Annual  Report.1  The  arrangement  used  the  1100-1300  nm  light 
from  a  Cnforsterite  laser  to  illuminate  the  sample,  a  Fourier  space  gate, 6  and  a  polarization  gate7  to  sort  out 
photons  that  carry  direct  image  information,  and  an  InGaAs  NIR  area  camera  to  record  2-D  transillumination 
images.  Spectroscopic  imaging  was  an  important  new  development  in  the  project  that  was  not  in  the  original 
SOW,  and  it  enabled  us  to  explore  and  demonstrate  the  diagnostic  potential  of  the  optical  imaging  that  was 
contemplated  in  TO  6  and  outlined  in  Task  13.  Results  of  spectroscopic  imaging  were  presented  in  earlier 
Annual  Reports1'3  and  will  be  discussed  later  in  this  Final  Report. 

5.2.  Time-sliced  Direct  Imaging  of  ex  vivo  Normal  and  Cancerous  Human  Breast  Tissues 

The  efficacy  of  the  experimental  arrangements  were  tested  and  evaluated  through  imaging  experiments  on  ex 
vivo  human  breast  tissue  samples  (TO  5-7,  Tasks  12-17).  The  samples  were  obtained  from  our  collaborators  at 
the  Memorial  Sloan  Kettering  Cancer  Center  (MSKCC),  as  well  as  from  National  Disease  Research  Interchange 
(NDRI)  under  Internal  Review  Board  (IRB)  approval  at  CCNY  and  were  characterized  by  biopsy  (Task  12)  and 
later  by  pathology  (Task  17).  Breast  tissue  specimens  from  patients  of  different  ages  (between  40-70  years)  and 
different  types  of  breast  cancers  (infiltrating  ductal  carcinoma,  infiltrating  lobular  carcinoma,  mucinous 
carcinoma,  invasive  ductal  carcinoma  etc.)  were  investigated  along  with  normal  tissues  from  the  same  patients. 
The  thickness  of  breast  tissue  specimens  normally  varied  between  5-10  mm  for  direct  imaging. 

The  tissue  specimens  were  held  between  glass  plates  under  compression  to  maintain  uniform  thickness  and  a 
sequence  of  time-sliced  2-D  transillumination  images  were  recorded  for  different  positions  of  the  time  gate. 

The  compression  proved  beneficial  (Task  13)  for  maintaining  uniform  thickness,  and  avoiding  voids  and  gaps 
in  the  specimen.  Time-sliced  imaging  measurements  of  tissues  with  tumor  (Task  14)  demonstrated  that  2-D 
transillumination  images  recorded  with  early  arriving  light  highlighted  the  tumor,  while  later  arriving  light 
accentuated  the  normal  tissue.  This  behavior  was  observed  for  different  types  of  breast  cancers  studied  in  the 
project.  Two-dimensional  time-sliced  images  showing  the  difference  between  normal  and  cancerous  tissues 
were  presented  in  Annual  Reports  1-3,  and  appeared  in  published  articles8'11  and  presentations.  Details  of  the 
experimental  parameters,  procedures,  and  results  are  presented  in  those  primary  sources  [Reference  8  and 
Reference  9  are  Appendix  1  and  Appendix  2,  respectively  of  First  Annual  Report,1  Reference  10  is  Appendix  1 
of  Third  Annual  Report,2  and  Reference  11  is  Appendix  5  of  this  Final  Report]. 

5.3.  Spectroscopic  Direct  Imaging  of  ex  vivo  Normal  and  Cancerous  Human  Breast 
Tissues 

Spectroscopic  imaging  experiments  were  carried  out  on  the  ex  vivo  breast  tissue  samples  mentioned  in  Section 
5.3.  In  most  cases,  time-sliced  and  spectroscopic  imaging  measurements  were  carried  out  in  a  coordinated 
manner  so  that  the  results  from  the  two  approaches  could  be  compared.  Two-dimensional  transillumination 
images  were  recorded  using  light  of  several  wavelengths,  such  as,  1210, 1225, 1235, 1250, 1280,  and  1300  nm 
[TO  5-7,  Tasks  13  and  15]  from  the  Cnforsterite  laser.  As  a  proof  of  principle  of  the  diagnostic  ability  of 
spectroscopic  imaging,  we  first  carried  out  measurements  on  normal  ex  vivo  breast  tissues  with  adipose  and 
fibrous  regions.4,5’12  Transillumination  images  were  recorded  with  light  whose  wavelength  is  close  to  the 
adipose  absorption  resonance13  at  1203  nm  and  compared  with  those  recorded  with  off-resonance  wavelength. 
Images  recorded  with  near-resonant  light  distinguished  the  adipose  region  with  higher  contrast  indicating  the 
diagnostic  potential  of  the  optical  imaging,  if  appropriate  wavelengths  could  be  identified. 

While  such  a  marked  contrast  was  not  observed  between  the  normal  and  cancerous  regions  in  the  2-D 
transillumination  images  of  the  breast  tissue  specimens,  two  salient  features  emerged.  There  was  an  overall 


5 


*  higher  light  transmisision  through  the  cancerous  region  at  all  of  the  above-mentioned  wavelengths.  What  is 
even  more  interesting,  a  wavelength-dependent  difference  in  light  transmission  through  the  cancerous  and 
normal  tissues  was  observed.10,11  The  tumor-to-normal  intensity  ratio,  /?TN(^)=ftumor(^)/4ormai(^)5  where, 

/tun,or(X)  and  /normal),  are  averaged  transmitted  intensities  through  the  tumor  and  the  normal  tissue  regions, 
respectively,  had  a  value  of  2.2  at  1225  nm,  and  1.5  at  1300  nm,  a  significant  difference.11 

In  an  alternative  application  involving  Warthin’s  tumor  of  the  parotid  gland,  the  value  of  was  found  to  be  2.9  at 
1225  nm  and  1.7  at  1285  nm.14  Implications  of  these  wavelength-dependent  transmission  ratios  are  further 
discussed  in  these  appended  references  [Reference  1 1  and  Reference  14  are  Appendix  5  and  Appendix  2, 
respectively,  of  this  Final  Report]. 

5.4  Development  of  Analytical  Solution  of  Radiative  Transfer  Equation 

Development  of  inverse  image  reconstruction  methods  requires  a  theoretical  model  that  provides  an  accurate 
description  of  photon  transport  through  highly  scattering  turbid  media.  The  commonly  used  forward  model  for 
light  propagation  in  breast  tissue  is  based  on  the  Diffusion  Approximation15  (DA)  of  the  Radiative  Transfer 
Equation  (RTE).  DA-based  forward  models  cannot  account  for  ballistic  and  snake  photons,  and  lead  to  large 
errors  in  the  region  near  sources,  where  the  weight  function  is  most  important.  We  have  made  a  major 
breakthrough  in  developing  exact  analytical  solutions  of  the  RTE  and  went  on  to  build  forward  models  based  on 
that  cumulant  solution,  as  a  part  of  our  endeavor  to  develop  novel  inverse  image  reconstruction  methods  (TO  4, 
Task  10;  TO  8,  Task  18).  Our  accomplishments  include  the  following. 

■  Development  of  an  analytical  solution  of  RTE,  based  on  cumulant  expansion,  in  an  infinite  uniform 
medium  with  an  arbitrary  phase  function.16'17,18  It  provides  an  explicit  analytical  expression  for  photon 
density  distribution  function  /( r,  s,  t),  as  a  function  of  position  r,  direction  of  light  s,  and  time  t.  The 
mean  position  and  the  half-width  at  half  maximum  height  (HWHM)  of  the  distribution  are  always  exact. 
Detailed  description  of  the  approach  was  presented  in  the  Second  Annual  Report,2  and  appended 
publication  [Appendix  2  and  Appendix  3  of  the  Second  Annual  Report2]. 

■  Extension  of  analytical  cumulant  solution  to  include  the  light  polarization,  i.e.,  solution  of  polarized 
photon  transport  equation  in  an  infinite  uniform  medium  that  can  take  into  account  the  vector  nature  of 
electromagnetic  radiation  and  can  handle  the  polarization  property  of  light.19  This  highly  significant 
advance  was  presented  in  the  Third  Annual  Report,3  and  publication  was  appended  [Appendix  2  in  Third 
Annual  Report3]. 

■  Development  of  the  linear  cumulant  forward  model  (CFM)  based  on  cumulant  analytical  solution  of  the 
RTE.20"21  This  CFM  may  be  used  with  time-resolved,  continuous-wave,  and  frequency-domain  data  to 
reconstruct  images  that  are  much  more  accurate  than  that  provided  by  DA-based  models.  [Reference  20 
and  Reference  21  are  Appendix  1  and  Appendix  4,  respectively,  of  this  Final  Report.] 

5.5  Development  and  Demonstration  of  Inverse  Image  Reconstruction 

We  followed  several  schemes  to  develop  and  implement  inverse  image  reconstruction  (HR)  methods, (TO  4, 
Tasks  8, 9)  and  demonstrate  have  developed  algorithms  for  realizing  fast  three-dimensional  inverse  image 
reconstruction  (TO  4,  Tasks  10, 11;  and  TO  8,  Task  18).  Included  in  the  development  and  demonstration  of 
HR  methods  are  the  following. 

■  Development  of  an  DR  formalism  that  used  a  point  source  illumination,  a  sequence  of  2-D  time-sliced 
images  recorded  with  a  time-gated  CCD  as  input  data,  a  DA-based  forward  model  that  assumes  a  slab 
geometry  in  cylindrical  coordinates,  and  a  Green’s  function  perturbative  approach  under  the  Rytov 
approximation  for  obtaining  a  linear  inversion  algorithm.22  A  salient  feature  of  the  formalism  is  the  use 
of  a  2-D  matrix  inversion  with  a  1-D  Fourier  transform  inversion  based  on  the  cylindrical  symmetry  that 
greatly  reduced  the  computation  time.  [Detailed  in  Appendix  3  of  the  First  Annual  Report1] 

■  Development  of  an  DR  approach  based  on  the  concept  of  propagation  of  spatial  Fourier  components  of 
the  scattered  wave  field  that  can  use  both  transmitted  and  backscattered  light  from  the  scattering 
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medium  containing  absorbing  and  scattering  inhomogeneities  to  construct  3-D  images.23  [Detailed  in 
Appendix  4  of  the  Second  Annual  Report2  and  Appendix  4  of  the  Third  Annual  Report3] 

*  Demonstration  of  3-D  image  reconstruction  based  on  the  diffusion  forward  model.22,23  Further 
improvement  of  these  inverse  algorithms  by  combining  with  the  forward  model  based  on  solution  of 
RTE  is  being  pursued.20,21  Details  of  this  development  appear  in  References  20  and  21  which  are 
included  as  Appendix  1  and  Appendix  4  of  this  Final  Report,  respectively. 

5.6.  Correlation  with  Existing  Methods 

It  is  imperative  that  any  new  detection  and  diagnostic  breast  imaging  technology  be  validated  by  comparison 
with  existing  methods.  Biopsy,  x-ray,  nuclear  magnetic  resonance  (NMR),  and  pathology  are  among  the 
commonly  used  methods  for  tumor  detection  and  diagnosis.  We  have  chosen  breast  tissue  specimens  for  study 
based  on  biopsy  results,  and  following  optical  imaging  obtained  histopathological  evaluation  of  the  samples 
(TO  5-7,  Task  17)  for  understanding  and  validating  results  of  the  optical  approach.  Comparison  with 
histopathology  provided  insight  into  the  composition  and  nature  of  tumor  and  normal  tissues,  and  helped 
explain  observed  optical  results  (See,  for  example,  Reference  14  which  is  Appendix  2  of  this  Final  Report). 

While  biopsy  and  pathology  provide  useful  information  for  ex  vivo  specimens,  alternative  approaches  are 
needed  for  validation  of  in  vivo  studies.  NMR  can  provide  an  in  vivo  assessment.  So,  we  started  measurements 
on  the  same  excised  breast  tissue  specimens  using  both  optical  and  NMR  techniques  to  investigate  how  the 
results  correlate  (TO  5-7,  Task  17)  and  reported  the  initial  results  in  the  Third  Annual  Report.  We  continued 
comparison  with  NMR  measurements  carried  out  at  the  Memorial  Sloan  Kettering  Cancer  Center  by  our 
collaborator  Dr.  Jason  Koutcher.  NMR  Ti-weighted  and  T2-weighted  images  of  the  breast  tissue  specimens 
exhibited  good  correlation  with  optical  measurements.  The  results  were  presented  in  a  major  conference,  and  a 
manuscript  has  been  submitted  for  publication  in  the  conference  proceedings  (Appendix  5  of  this  Final  Report). 
We  refer  to  the  appended  manuscript  for  further  details. 

Comparison  of  the  results  from  direct  imaging  and  image  reconstruction  (TO  5-7,  Task  16)  was  limited  by  the 
fact  that  beyond  a  tissue  thickness  of  10  mm  direct  imaging  did  not  produce  well  resolved  images,  while  DA- 
based  reconstruction  required  thicker  tissue.  However,  the  differences  observed  in  the  direct  time-sliced  and 
spectroscopic  images  provided  the  basis  for  the  reconstruction  images  to  be  viable. 

(6)  KEY  RESEARCH  ACCOMPLISHMENTS 

•  Demonstrated  that  time-sliced  2-D  transillumination  images  recorded  with  earlier  temporal  slices  of 
transmitted  light  highlight  tumor  while  those  recorded  with  later  slices  accentuate  normal  fibrous 
tissues  in  ex  vivo  human  breast  tissue  specimens. 

•  Near-infrared  2-D  spectroscopic  imaging  experiments  on  more  breast  tissue  specimens  could 
distinguish  between  normal  adipose  and  fibrous  tissues  when  transillumination  images  were  recorded 
with  light  near  resonant  with  the  adipose  absorption  at  1203  nm  and  compared  with  that  recorded  with 
off-resonance  light.  Extending  the  measurements  to  normal  and  cancerous  breast  tissues  we  observe 
that  the  ratio  R  of  light  intensity  transmitted  through  the  cancerous  tissue  to  that  through  the 
corresponding  normal  tissue  as  a  function  of  wavelength  might  be  used  as  a  useful  parameter  for 
cancer  identification. 

•  Developed  analytical  solutions  of  the  scalar  and  vector  forms  of  the  radiative  transport  equation  based 
on  a  cumulant  expansion  that  enable  much  more  accurate  description  of  photon  transport  thorough 
scattering  media  than  that  afforded  by  the  Diffusion  Approximation.  Cumulant  solutions  can  describe 
all  scattered  photons  including  snake  photons,  and  the  solutions  of  the  vector  RTE  can  take  into 
consideration  the  polarization  properties  of  light  transmitting  through  a  highly  scattering  medium. 
Extended  Cumulant  Solution  of  scalar  transport  equation  to  semi-infinite  and  slab  geometries  making 
it  more  suited  for  practical  application  that  involves  finite  boundaries.  These  analytical  solutions 
provide  basis  for  developing  more  accurate  forward  models  and  inverse  image  reconstruction 
algorithms. 
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•  Developed  an  inversion  algorithm  that  uses  a  temporal  sequence  time-sliced  2-D  images  as  input  data, 
a  DA-based  forward  model  and  a  perturbative  linear  inversion  approach,  and  demonstrated  fast  3-D 
image  reconstruction. 

•  Developed  a  forward  model  based  on  the  cumulant  solution  of  the  RTE  for  use  in  more  sophisticated 
inversion  algorithms  for  better  image  reconstruction  of  tumors  in  breast,  and  objects  in  obscuring 
media,  in  general. 

•  Demonstrated  correlation  between  results  of  time-sliced  optical  imaging  and  NMR  imaging  methods. 

(7)  REPORTABLE  OUTCOMES 

7.1.  Articles 

1.  S.  K.  Gayen  and  R.  R.  Alfano,  “Sensing  lesions  in  tissues  with  light,”  Opt.  Express  4, 475  (1999). 

2.  W.  Cai,  S.  K.  Gayen,  M.  Xu,  M.  Zevallos,  M.  Alrubaiee,  M.  Lax,  and  R.  R.  Alfano,  "Optical  three- 
dimensional  inverse  image  reconstruction  of  objects  in  turbid  media  from  ultrafast  time-sliced  optical 
transmission  measurements,"  Appl.  Opt.  38, 4237  (1999). 

3.  S.  K.  Gayen,  M.  E.  Zevallos,  B.  B.  Das,  and  R.  R.  Alfano,  "Near-infrared  spectroscopic  and  time-sliced 
imaging  of  human  breast  tissues,"  SPIE  Proceedings  on  Optical  Tomography  and  Spectroscopy  of  Tissue 
fflVol.  3597, 508  (1999). 

4.  M.  Xu,  S.  K.  Gayen,  W.  Cai,  M.  E.  Zevallos,  M.  Lax,  and  R.  R.  Alfano,  "Time-sliced  three-dimensional 
inverse  image  reconstruction  of  objects  in  highly  scattering  media."  SPIE  Proceedings  on  Optical 
Tomography  and  Spectroscopy  of  Tissue  HI  Vol.  3597, 2  (1999). 

5.  M.  Zevallos,  S.  K.  Gayen,  B.  B.  Das,  M.  Alrubaiee,  R.  R.  Alfano,  “Imaging  of  bones  inside  tissues  using  a 
picosecond  electronic  time  gate,”  IEEE  J.  Select.  Topics  Quantum  Electron.,  5, 916  (1999). 

6.  W.  Cai,  M.  Lax,  and  R.  R.  Alfano,  “Analytical  solution  of  the  elastic  Boltzmann  transport  equation  in  an 
infinite  uniform  medium  using  cumulant  expansion,”  J.  Phys.  Chem.  B  104, 3996  (2000). 

7.  W.  Cai,  M.  Lax,  and  R.  R.  Alfano,  “Cumulant  solution  of  the  elastic  Boltzmann  transport  equation  in  an 
infinite  uniform  medium,”  Phys.  Rev.  E  61, 3871  (2000). 

8.  S.  K.  Gayen,  M.  Alrubaiee,  M.  E.  Zevallos  and  R.  R.  Alfano,  ‘Temporally  and  spectrally  resolved  optical 
imaging  of  normal  and  cancerous  human  breast  tissues,”  in  the  Proceedings  of  the  Inter-Institute 
Workshop  on  In  Vivo  Optical  Imaging  at  the  NIH,  Amir  Gandjbakhche,  ed.  (Optical  Society  of  America, 
Washington,  DC,  2000),  pp.  142-147. 

9.  M.  Xu,  M.  Lax,  and  R.  R.  Alfano,  ‘Time-resolved  optical  diffuse  tomography,”  J.  Opt.  Soc.  Am.  A  18, 
1535(2001). 

10.  W.  Cai,  M.  Lax,  R.  R.  Alfano,  “Analytical  solution  of  the  polarized  photon  transport  equation  in  an 
infinite  uniform  medium  using  cumulant  expansion,”  Phys.  Rev.  E  63, 016606-1(2001). 

11.  M.  Xu,  W.  Cai,  M.  Lax,  R.  R.  Alfano,  “Photon  transport  forward  model  for  image  in  turbid  media”  Opt. 
Lett.  26  1066-1068  (2001). 

12.  S.  K.  Gayen,  M.  Alrubaiee,  H.  E.  Savage,  S.  P.  Schantz,  and  R.  R.  Alfano,  “Parotid  gland  tissues 
investigated  by  picosecond  time-gated  and  spectroscopic  imaging  techniques,”  J.  Select.  Topics  Quantum 
Electron.  7, 906  (2001). 

13.  W.  Cai,  M.  Xu,  M.  Lax,  R.  R.  Alfano,  “Diffusion  coefficient  depends  on  time  not  on  absorption”  Opt. 
Lett.  27  731-733  (2002). 

14.  M.  Xu,  W.  Cai,  M.  Lax,  R.  R.  Alfano,  “Photon  migration  in  turbid  media  using  a  cumulant  approximation 
to  radiative  transfer”  Phys.  Rev.  E  65  066609  (2002). 

15.  M.  Al-Rubaiee,  S.  K.  Gayen,  J.  A.  Koutcher,  and  R.  R.  Alfano,  "Cancerous  and  normal  human  tissues 
investigated  by  near-infrared  time-resolved  and  spectroscopic  imaging  techniques."  SPIE  Proceedings  on 
Optical  Tomography  and  Spectroscopy  of  Tissue  (to  be  published). 
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7.2.  Abstracts  and  Presentations 

1.  M.  Xu,  S.  K.  Gayen,  W.  Cai,  M.  E.  Zevallos,  M.  Lax,  and  R.  R.  Alfano,  "Time-sliced  three-dimensional 
inverse  image  reconstruction  of  objects  in  highly  scattering  media."  Paper  3597-02  presented  at  the  SPIE's 
International  Symposium  on  Biomedical  Optics,  BiOS  '99/  Photonics  West,  23-29  January,  San  Jose, 
California. 

2.  S.  K.  Gayen,  M.  E.  Zevallos,  B.  B.  Das,  and  R.  R.  Alfano,  "Near-infrared  spectroscopic  and  time-sliced 
imaging  of  human  breast  tissues."  Paper  3597-84  presented  at  the  SPIE's  International  Symposium  on 
Biomedical  Optics,  BiOS  '99/  Photonics  West,  23-29  January,  San  Jose,  California. 

3.  S.  K.  Gayen,  M.  E.  Zevallos,  and  R.  R.  Alfano,  "Time-resolved  transillumination  imaging  of  normal  and 
cancerous  human  breast  tissues  with  a  picosecond  electronic  time  gate,"  Bull.  Am.  Phys.  Soc.  44,  117 
(1999).  Paper  BC32  7  presented  at  the  American  Physical  Society  Centennial  Meeting,  March  20-26, 
1999,  Atlanta,  Georgia. 

4.  S.  K.  Gayen,  M.  E.  Zevallos,  and  R.  R.  Alfano,  "Time-sliced  and  spectroscopic  correlated  two- 
dimensional  near-infrared  imaging  for  cancer  detection."  Paper  CFL4  presented  at  the  Conference  on 
Lasers  and  Electro-Optics  (CLEO’99),  May  23-28, 1999,  Baltimore,  Maryland. 

5.  S.  K.  Gayen,  M.  Alrubaiee,  R.  R.  Alfano,  J.  Koutcher,  and  H.  Savage,  "Electronic  time-gated  and 
spectroscopic  near-infrared  imaging  of  lesions  in  human  tissues,"  Bull.  Am.  Phys.  Soc.  45,  958  (2000). 
Paper  Y13  10  presented  at  the  March  Meeting  of  the  American  Physical  Society,  20-24  March  2000, 
Minneapolis,  MN. 

6.  S.  K.  Gayen,  M.  Alrubaiee,  and  R.  R.  Alfano,  ‘Time-sliced  and  spectroscopic  two-dimensional  imaging 
of  normal  and  cancerous  human  breast  tissues,”  The  Department  of  Defense  Breast  Cancer  Research 
Program  Meeting  Era  of  Hope,  Proceeding  Volume  I,  p.195.  Presented  at  the  Hilton  Atlanta  and  Towers, 
Atlanta,  GA,  June  8-11, 2000. 

7.  M.  Xu,  M.  Lax,  and  R.  R.  Alfano,  "Time-resolved  Fourier  diffuse  optical  tomography,"  in  the  Technical 
Digest  of  Biomedical  Topical  Meetings,  TOPS  Volume  38,  Optical  Society  of  America  (2000)  pp.  345- 
347.  Paper  TuC4  presented  at  the  Biomedical  Topical  Meetings,  Miami  Beach,  Florida,  April  2-5, 2000. 

8.  W.  Cai,  M.  Lax,  and  R.  R.  Alfano,  "Analytical  solution  of  the  polarized  photon  transport  equation  in  an 
infinite  uniform  scattering  medium,"  in  the  Technical  Digest  of  Biomedical  Topical  Meetings,  TOPS 
Volume  38,  Optical  Society  of  America  (2000)  pp.  62-64.  Paper  TuC4  presented  at  the  Biomedical 
Topical  Meetings,  Miami  Beach,  Florida,  April  2-5, 2000. 

9.  M.  Alrubaiee,  S.  K.  Gayen,  H.  E.  Savage,  S.  P.  Schantz,  and  R.  R.  Alfano,  "Picosecond  time-gated  and 
optical  spectroscopic  imaging  of  parotid  gland  tissues."  Paper  FP-2  presented  at  Frontiers  of  Photonics 
Symposium,  City  College  of  New  York,  November  5, 2001,  New  York. 

10.  M.  Xu,  W.  Cai,  M.  Lax,  R.  R.  Alfano,  “Prior  information  and  noise  in  three-dimensional  optical  image 
reconstruction,”  in  the  Technical  Digest  of  OSA  Biomedical  Topical  Meetings  (2002)  p30.  Paper  SuB7-l 
presented  at  the  Biomedical  Topical  Meetings,  Miami  Beach,  Florida,  April  7-10,  (2002). 

11.  M.  Xu,  W.  Cai,  R.  R.  Alfano,  ‘Three-dimensional  hybrid-dual-Fourier  tomography  in  turbid  media  using 
multiple  sources  and  multiple  detectors,”  Paper  presented  at  the  Third  Inter-Institute  Workshop  on 
Diagnostic  Optical  Imaging  and  Spectroscopy:  The  Clinical  Adventure  on  NIH  campus  in  Bethesda,  MD, 
Sept.  26-27  (2002). 

12.  W.  Cai,  S.  K.  Gayen,  M.  Xu,  R.  R.  Alfano,  “Improving  inverse  reconstruction  problem  for  three- 
dimensional  optical  image  of  breast,”  The  Department  of  Defense  Breast  Cancer  Research  Program 
Meeting  Era  of  Hope,  proceedings  Vol  HI  p48-l,  presented  at  Orange  County  Convention  Center, 
Orlando,  Florida,  Sept  25-28  (2002). 

13.  M.  Al-Rubaiee,  S.  K.  Gayen,  W.  Cai,  M.  Xu,  J.  A.  Koutcher,  and  R.  R.  Alfano,  "Near-infrared  photonics 
imaging  of  human  breast  tissue."  Paper  P48-2  presented  at  Era  of  Hope,  Department  of  Defense  Breast 
Cancer  Research  Program  Meeting,  September  25-28, 2002,  Orlando,  Florida. 

14.  M.  Al-Rubaiee,  S.  K.  Gayen,  J.  A.  Koutcher,  and  R.  R.  Alfano,  "Cancerous  and  normal  human  tissues 
investigated  by  near-infrared  time-resolved  and  spectroscopic  imaging  techniques."  Paper  4955-28 
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presented  at  the  Optical  Tomography  and  Spectroscopy  of  Tissue  V  Conference  of  SPIE’s  Photonics  West 
BiOS  2003/  Photonics  West,  January  25-31, 2003,  San  Jose,  California. 

7.3.  Patents 

1.  R.  R.  Alfano,  W.  Cai,  S.  K.  Gayen,  ‘Time-resolved  Diffusion  Tomographic  2D  and  3D  Imaging  in  Highly 
Scattering  Turbid  Media”  U.S.  Patent  6,108,576,  issued  August  22, 2000. 

2.  R.  R.  Alfano,  W.  Cai,  M.  Lax,  ‘Time-resolved  optical  backscattering  tomographic  image  reconstruction  in 
scattering  turbid  media”  U.  S.  Patent  6,205,353,  issued  March  20, 2001. 

7.4.  Degrees  Obtained 

1.  Manuel  Zevallos  received  his  Ph.  D.  in  Electrical  Engineering  in  1999.  Part  of  the  research  for  his  thesis 
entitled,  “Near-infrared  optical  imaging  and  light  propagation  in  highly  scattering  random  media,”  was 
supported  by  this  grant. 

2.  Min  Xu  received  his  Ph.  D.  in  Physics  in  2001.  He  was  supported  in  part  by  this  grant.  The  title  of  his 
thesis  is,  “Optical  image  reconstruction  in  highly  scattering  media.” 

7.5.  Funding  Applied  for  Based  on  Work  Supported  by  this  Award 

1.  ‘Two-dimensional  and  three-dimensional  optical  mammography,”  was  submitted  to  USAMRMC  on  6/1 1/02 
in  response  to  BCRP-2002. 

2.  ‘Time-resolved  2-D  and  3-D  imaging  of  breast  and  Prostate,”  submitted  to  the  Network  for  Translational 
Research:  Optical  Imaging  program  of  NCI  on  2/26/03. 

3.  “Center  for  Photonic  Imaging  and  Diagnosis  of  Breast  Cancer,”  a  pre-proposal  submitted  to  USAMRMC, 
DoD  2003  Breast  Cancer  Research  Program.  We  have  recently  been  notified  to  submit  a  full  proposal. 

7.6.  Research  Opportunities  Applied  for  and/or  Received  Based  on  Experience/Training 
Supported  by  this  Award 

Dr.  Min  Xu  was  supported  in  part  by  this  grant  during  his  doctoral  work.  Based  on  the  training  he  received,  he 
went  on  to  apply  for  a  Postdoctoral  Traineeship  Award.  His  proposal,  ‘Time-resolved  spectral  optical  breast 
tomography,”  was  supported  by  USAMRMC  (Grant:  DAMD 17-02- 1-05 16)  for  the  total  amount  of  $150,000 
covering  the  period  of  May  15, 2002-  May  14, 2005. 

7.7.  Personnel  Supported  by  the  Research  Effort 

Persons  Directly  Supported 

1.  R.  R.  Alfano:  As  the  Principal  Investigator  (PI)  of  the  project  R.  R.  Alfano,  Distinguished  Professor  of 
Science  and  Engineering,  was  involved  in  the  overall  supervision  of  the  research  carried  out. 

2.  W.  Cai,  Research  Associate:  Pursued  development  of  analytical  solutions  of  the  Radiative  Transfer 
Equation  and  algorithms  for  inverse  image  reconstruction  as  a  research  associate. 

3.  B.  B.  Das,  Research  Associate:  Supported  in  part  by  this  grant,  B.  B.  Das,  was  involved  in  the 
measurement  of  the  temporal  profiles  of  light  pulses  transmitted  through  and  backscattered  from  breast 
tissues  and  highly  scattering  model  media  to  estimate  light  transport  parameters,  such  as,  transport  mean 
free  path  and  scattering  length. 

4.  S.  K.  Gayen:  As  a  Co-PI  of  the  project,  S.  K.  Gayen  was  responsible  for  directing  the  experimental 
effort,  and  coordinating  with  the  subcontractor,  and  the  theoretical  group.  He  started  as  a  Research 
Associate  and  went  on  to  become  an  Associate  Professor  of  Physics  during  the  period  covered  by  the 
project.  He  was  supported  in  part  by  the  project. 

5.  X.  Liang:  Responsibilities  included  maintenance  of  the  Ti.sapphire  laser  system,  and  help  carry  out 
imaging  experiment.  Received  partial  support  as  a  research  technician. 
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6.  Min  Xu:  Min  Xu  was  involved  in  the  development  of  image  reconstruction  algorithms,  and  extension  of 
the  cumulant  solution  of  radiative  transport  equation.  He  was  supported  in  part  as  a  graduate  student. 

He  received  his  Ph.  D.  and  went  on  to  earn  a  postdoctoral  trainee  grant  from  USAMRMC. 

Persons  Supported  Indirectly 

7.  M.  Alrubaiee:  M.  Alrubaiee  is  pursuing  his  graduate  work  in  the  area  of  research  pursued  by  this 
project,  carried  out  many  of  the  project  experiments.  He  did  not  receive  any  stipend  from  this  effort,  but 
the  project  was  supported  by  the  grant. 

8.  M.  E.  Zevallos:  M.  E.  Zevallos  received  his  Ph.  D.  in  1999,  and  part  of  his  thesis  work  was  done  on  this 
project.  His  stipend  was  derived  from  another  source. 

(8)  CONCLUSIONS 

The  work  carried  out  on  this  project  provides  important  information  about  the  development  of  time-resolved 
optical  and  spectroscopic  imaging  techniques  for  detection  and  diagnosis  of  breast  cancer.  First,  direct  time- 
sliced  2-D  transillumination  images  can  sort  out  tumor  from  normal  breast  tissue.  Images  recorded  with  earlier 
temporal  slices  of  transmitted  light  highlighted  cancerous  tissues  while  those  recorded  with  later  slices 
accentuated  normal  fibrous  tissues.  Second,  spectroscopic  imaging  shows  tissue  selectivity,  and  diagnostic 
potential.  Wavelength  dependence  of  a  ratio,  R  of  light  intensity  transmitted  through  the  tumor  and  normal 
tissue  shows  the  potential  to  be  used  as  a  useful  parameter  for  cancer  identification.  Third,  a  sequence  of  time- 
sliced  2-D  images  provides  a  wealth  of  information  for  3-D  inverse  image  reconstruction,  that  are  not 
commonly  available  from  other  optical  methods.  Fourth,  analytical  solutions  of  the  scalar  and  vector  photon 
transport  equations  using  Cumulant  Expansion  enables  more  complete  description  of  diffusive,  snake  and 
ballistic  photons  than  the  currently  used  methods,  and  can  accommodate  for  polarization  characteristics  of  light 
emerging  from  a  scattering  medium  and  will  be  useful  for  developing  better  image  reconstruction  approaches. 
Fifth,  the  theoretical  formalism  and  computer  algorithm  for  inverse  image  reconstruction  scheme  using  both 
forward  and  backscattered  light  shows  (with  simulated  data)  the  potential  to  provide  fast  3-D  images  of 
scattering  and  absorbing  objects  at  various  depths  inside  a  scattering  medium.  Finally,  the  good  agreement 
between  results  of  time-sliced  optical  measurement  and  NMR  measurements  indicates  the  possibility  of  NMR 
co-registration  for  validation  of  in  vivo  optical  imaging  measurements. 

We  plan  to  build  on  these  developments  and  pursue  in  vivo  breast  imaging  on  volunteers  following  testing  of 
the  techniques  on  model  breast  structures  made  from  excised  tissues.  We  submitted  a  White  Paper  entitled, 
“Center  for  Photonic  Imaging  and  Diagnosis  of  Breast  Cancer,”  in  response  to  DoD  Fiscal  Year  2003  Breast 
Cancer  Research  Program.  We  have  recently  been  invited  to  submit  a  full  proposal.  The  research  programs  of 
the  envisioned  Center  will  address  the  overarching  question  concerning  early  detection  and  diagnosis  of  breast 
cancer,  and  will  be  a  natural  culmination  of  the  work  carried  our  under  this  grant  support. 

“So  What  Section” 

The  implication  of  the  first  conclusion  is  that  time-sliced  imaging  offers  the  possibility  of  highlighting 
cancerous  lesions  in  human  breast.  The  second  conclusion  indicates  the  diagnostic  potential  of  optical  imaging 
based  upon  multi-spectral  measurements.  X-ray  mammography,  most  often  used  method,  cannot  diagnose 
cancer.  The  third,  fourth  and  fifth  conclusions  together  present  the  possibility  of  developing  robust  3-D  inverse 
image  reconstruction  formalisms,  that  in  addition  to  being  applicable  for  optical  mammography,  will  be  useful 
for  imaging  tumors  in  other  body  organs,  such  as,  prostate  and  thyroid  glands,  as  well  as,  objects  inside  other 
types  of  scattering  media,  such  as,  cloud,  fog,  smoke,  and  murky  water.  The  final  conclusion  points  to  the 
complementary  nature  of  optical  and  NMR  imaging,  and  the  possibility  of  using  them  together  for  better 
specificity.  Overall,  the  time-sliced  and  spectroscopic  imaging  approach  developed  and  pursued  in  this  project 
may  lead  to  the  “holy  grail”  for  obtaining  breast  cancer  detection  and  diagnostic  information  based  on 
wavelength  selection  and  time-sliced  image  reconstruction. 
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A  photon-transport  forward  model  for  image  reconstruction  in  turbid  media  is  derived  that  treats  weak  in¬ 
homogeneities  through  a  Born  approximation  of  the  Boltzmann  radiative  transfer  equation.  This  model  can 
conveniently  replace  the  commonly  used  diffusion  approximation  in  optical  tomography.  An  analytical  ex¬ 
pression  of  the  background  Green’s  function  is  obtained  from  the  cumulant  solution  of  the  Boltzmann  equation. 

Our  model  provides  the  correct  behavior  of  photon  migration  at  early  times  and  reduces  at  long  times  to  the 
center-moved  diffusion  approximation.  Numerical  comparisons  between  this  model  and  the  standard  and 
center-moved  diffusion  models  are  presented.  ©  2001  Optical  Society  of  America 
OCIS  codes:  170.6960,  290.7050,  030.5620. 


Probing  the  internal  properties  of  highly  scattering 
turbid  media  with  photons  has  a  variety  of  applications 
in  geophysics,  radio  astronomy,  and  medical  tomogra¬ 
phy.  Imaging  based  on  the  diffusion  approximation 
has  been  pursued  over  the  past  decade  because  it  cap¬ 
tures  the  core  characteristic  of  light  migration  in  tur¬ 
bid  media  and  is  easy  to  implement.1-3  However,  in 
the  diffusion  approximation,  light  is  assumed  to  dif¬ 
fuse  from  a  fixed  source  with  a  constant  diffusion  co¬ 
efficient  throughout  the  full  time  range  when  photons 
propagate  inside  the  uniform  medium.4  This  assump¬ 
tion  is  invalid  when  the  incident  photon  retains  its 
early  time  directionality  preference.  To  account  for 
this  difficulty,  a  common  practice  is  to  assume  that 
all  the  incident  photons  are  initially  scattered  at  a 
depth  zo  =  It  (transport  mean  free  path)  inside  the  tur¬ 
bid  medium,5  which  we  call  the  center-moved  diffusion 
model  (CDM).  But  CDM  breaks  the  reciprocity  theo¬ 
rem  and  still  fails  in  the  description  of  photon  propa¬ 
gation  at  early  times.6 

To  fully  account  for  photon  migration  in  a  turbid 
medium,  one  must  use  the  radiative  transfer  equation 
instead: 

-^-7(r,s,£)  +  cs  •  Vr7( r,  s,t)  +  c[>s(r)  +  /*a(r)]/(r,s,  t) 


=  cfjis( r)  j  ds'P(s,  s')7(r,  s',  f)ds'  +  q( r,  s,  *) ,  (1) 

where  7(r,  s,  t)  is  the  photon-distribution  function  de¬ 
pending  on  position  r,  direction  s,  and  time  t\  c  is  light 
speed  inside  the  medium;  fxa  and  fis  denote  the  posi¬ 
tion-dependent  absorption  and  scattering  coefficients; 
q( r,  s,  t)  is  the  photon-source  strength;  and  P(s,  s')  is 
the  normalized  phase  function  of  light  propagation  in 
the  medium. 

Recently,  we  derived  an  analytical  solution  for  the 
photon  distribution,  I^\r,  s,  t ),  and  photon  density, 
AT(0)( r,  t),  in  an  infinite  uniform  medium,  with  exact 
spatial  center  position  and  exact  spatial  half-width,  at 
any  direction  and  time7;  the  exact  solution  up  to  an 
arbitrary  order  of  cumulants  was  also  derived.8  The 
photon-density  distribution  is  found  to  have  a  center 
that  advances  in  time  and  an  ellipsoidal  contour  that 
grows  and  changes  shape,  providing  a  clear  picture 


of  the  time  evolution  of  light  migration  from  the  ini¬ 
tial  ballistic  to  the  final  diffusive  regime.  A  forward 
model  adding  inhomogeneity  to  the  analytical  expres¬ 
sion  of/(0)(r,  s,  t)  involves  a  complicated  numerical  in¬ 
tegration  over  angular  parameters.  In  this  Letter  we 
use  an  approximate  expression  of  7(0)(r,  s,  t)  that  not 
only  retains  the  main  features  of  photon  propagation  in 
both  early  and  long  time  limits  but  also  is  convenient 
for  building  a  forward  model  to  account  for  weak  in¬ 
homogeneities  of  the  medium  that  are  treatable  in  the 
Born  approximation. 

The  photon  distribution  in  an  infinite  uniform 
medium,  7(0)(r,  s,  t),  is  assumed  to  have  the  following 
form: 


7(0)(r,  s,  t |r0,  s0,  t0)  =  N(0\ r,  t |r0,  s0,  t0)F{0)(s,  t |s0,  to) 

-  f-D(t-  *o)8  •  Vr2V(0)(r,  f|ro,  s0,  to) ,  (2) 

4tt 

where  7V(0)(r,  f|r0,  s0,  fo)  is  the  photon  density  for 
a  point  pulse  propagating  along  so  at  position 
ro  and  time  to  in  an  infinite  uniform  medium, 
F(0)(s,  £|so,*o)  is  the  known  exact  photon  distri¬ 
bution  in  light-direction  space,  and  D(t  -  to) 
is  the  time-dependent  diffusion  coefficient.  The 
full  definitions  of  these  quantities  originated  in 
Ref.  7  and  are  given  as  follows:  F(0)(s,  t\&o,  to)  = 
(47 +  l)exp [~gi(t  -  t0)]Pi( s  •  s0),  where 
gi  =  cpts[  1  -  ai/(2l  +  1)],  ai  are  the  coefficients 
in  the  Legendre  expansion  of  the  phase  function, 
P(s,  s')  =  (477)"1Z/a/^(s  ■  s'),  especially,  go  =  0 
and  g\  =  cpLs  ,  where  fxs'  is  the  reduced  scattering 
coefficient.  The  diffusion  coefficient  D{t)  is  taken  to 
be  an  average  of  the  time-dependent  semimajor  axis, 
Dzz ,  and  semiminor  axes  Dxx  =  Dyy  of  the  diffusion 
coefficient  ellipsoid: 


Dxx  "f*  Dvv  D zz  c 

”<*>-  3  "at 


81 


1  -  exp(-gif) 
gi 2 


_  [1  -  exp(-gif)]2 

2gi2 

The  photon  density  is  given  by 


(3) 
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JV(0)(r,f|r0,so,<0) 


_ 1 _ 

[4  TrD(t  -  t0)c(t  -  f0)F2 


x  exp 


[r  -  r0  -  s0A(?  -  <0)]2 
4 D(t  -  to)c(t  -  £0) 


exp {~nac(t  -  f0)], 


(4) 

where  A (t)  =  c[l  -  exp(-gi£)]/gi  is  the  average  cen¬ 
ter  of  photons,  which  moves  with  speed  c  initially  and 
stops  at  c/gi  =  lt  in  the  long  time  limit. 

At  the  early  times  t  — *  to,  the  first  term  of  Eq.  (2) 
dominates,  and  F<0)(s,  <|s0,  f0)  —  S(s  -  s0),  D(t  -  t0)  -+ 
c2(t  -  tofgs'/d  —  0,  and  2V(0)(r,  t|r0,  s0,  to)  — 
5[r  -  r0  -  c(t  -  *o)so];  thus  /(0)( r,  s,  f|r0,s0,  t0 )  pro¬ 
vides  a  correct  picture  of  ballistic  motion  at  the  speed 
of  light  along  the  incident  direction,  s0.  In  the  long 
time  limit,  F(0)(s,  t\s0,  t0)  —  (47r)_1,  D(t )  —  (3/*s')~\ 
Eq.  (2)  reduces  to  the  photon  distribution  of  the  CDM 
approximation. 

A  perturbative  method  is  then  used  to  obtain  the 
forward  model  when  weak  inhomogeneities  are  intro¬ 
duced  into  the  otherwise  uniform  medium.  Making  a 
perturbation  expansion  of  the  radiative  transfer  equa¬ 
tion  (1)  to  the  first-order  Born  approximation,  we  de¬ 
rive  the  change  in  photon  distribution  from  7(0)(r,  s,  t), 
which  is  the  photon  distribution  of  the  uniform  back¬ 
ground,  as 


SJ(r,s,f|ro,s0,to) 

=  j  dt'  j  dr'  J  ds' I <0)(r',  -s',t  -  t'\r, -s) 
x  j J  ds"c8[fisP] (s',  s",  r')7(0)(r',  s",  t  -  «0|r0,s0) 


-  c[S/*s(r')  +  5/LtQ(r')]/(0)(r',s',<'  -  <0|r0,s0)  > 


(5) 


where  Sfia,  8fis,  and  8[/isP]  are  the  changes  of  the 
absorption  coefficient,  the  scattering  coefficient,  and 
the  angular-dependent  differential  scattering  coeffi¬ 
cient,  respectively,  from  the  background  to  the  inhomo¬ 
geneity.  The  optical  reciprocity  theorem,  7(0)(r,  s,  t  - 
t' Ir',  s')  =  7<0)(r',  -s',  t  -  t' |r,  -s),  is  used  to  obtain 
Eq.  (5). 

Expanding  <5[/xsP]  in  Legendre  polynomials  and 
substituting  Eq.  (2)  into  Eq.  (5),  the  integrations 
over  angular  variables  in  Eq.  (5)  can  be  analytically 
performed.  We  obtain 


5J(r,s,*|r0,So,*o)  =  J  dt'  J  dr'c8fia(r') 
X  N(0)(r',t  -  t'\r,-s)N{0\r',t  -  t0\r0,s0) 


dr 'D(t  -  t')D(t‘  - 10)  [S flair')  +Sfis  '(r')] 


X  Vr'Nm(r',t  -  t'\r,-S)-Vr,Nm(r',t'- t0 |r0,s0) 

+  4n  f  J  -<')[^Mo(r')+5/zs'(r')]exp(-^1<') 


X  {N(0)(r',t 


^olr,  ~s)s  •  Vr-A/'(0)(r',  t  -  P|r0,so) 


-  s0  •  Vr'Ni0)( r',  t  -  «'|r,  -s)N(0)(r',  t'  -  *0|r0,  s0»  (6) 


after  neglecting  fast-decaying  terms  involving 
exp(-2g;£)  for  l  >  1. 

The  photon-transport  forward  model  (PTFM), 
Eq.  (6),  is  the  main  result  in  this  Letter  and  obeys 
the  reciprocal  relation  8I( r,  s,  t  -  f0|r0,  s0)  = 
<5/(r0,  —So ,t  —  <0|r,  —  s).  In  the  long  time  limit,  the 
term  in  Eq.  (6)  that  contains  the  exponential  decay 
factor  exp(-gi£')  can  be  neglected,  and  the  change  in 
photon  density,  4n8I(r,  s,  t)  in  the  PTFM,  is  reduced 
to  the  diffuse  imaging  model  [Eq.  (14)  in  Ref.  9], 

To  show  the  difference  between  our  model  and 
the  diffusion  models,  we  consider  a  point  photon 
pulse  8(r)8(s  -  z)8(t)  propagating  inside  an  infinite 
scattering  turbid  medium  with  absorption  coefficient 
Ha  =  0,  reduced  scattering  coefficient  fis'  =  1  mm-1, 
anisotropy  g  =  0.9,  and  refractive  index  n  =  1.33. 
These  optical  properties  are  similar  to  those  of  a 
typical  breast  tissue.  A  Henyey-Greenstein  phase 
function  is  adopted  in  the  following  calculations.10 

The  time-resolved  profiles  of  transmission 
7(0)(r,  z,  t)  at  position  r  =  (0,  0,  5)  mm  and  backscat- 
tering  7(0)(r,  ~z,  t)  at  position  r  =  (0,  2,  0)  mm  are 
shown  by  the  solid  curves  in  Figs.  1(a)  and  1(b), 
respectively.  For  comparison,  the  photon  density 
divided  by  477  from  the  original  diffusion  model  (ODM) 
and  the  CDM  are  also  plotted. 

In  a  transmission  case  [Fig.  1(a)],  the  peak  photon 
intensity  in  the  ODM  arrives  too  late  compared  with 
the  experimental  result  reported  by  Yoo  et  al.,n 
whereas  in  the  CDM  the  artificial  adjustment  of  the 
source  position  for  lt  leads  to  an  arrival  of  photons 
faster  than  light  speed.  The  intensity  at  forward 
directions  from  our  model  is  stronger  than  that 
of  diffusion  models,  which  indicates  that  a  ceratin 
anisotropic  angular  distribution  remains  even  at  a 
distance  of  5 lt  from  the  source. 

In  the  case  of  backscattering  [Fig.  1(b)],  photons  dif¬ 
fuse  from  the  origin  (0,  0,  0)  to  (0,  2,  0)  mm  with  a 
constant  diffusion  coefficient  D  =  lt/ 3  in  the  ODM, 
whereas  photons  diffuse  from  the  adjusted  source  po¬ 
sition  (0,  0,  1)  mm  to  (0,  2,  0)  mm  with  the  same  con¬ 
stant  diffusion  coefficient  in  the  CDM.  The  photons 
in  our  model  are  backscattered  to  (0,  2,  0)  mm  later 
than  those  in  the  ODM  and  the  CDM  because  the 
center  of  photons  moves  forward  along  the  positive  z 
direction  and  diffuses  from  the  moving  center  with  a 
gradually  increasing  diffusion  coefficient  from  0  to  lt/3 
in  the  PTFM. 

Consider  a  forward  model  with  a  scattering  inho¬ 
mogeneity  8fis'  —  0.1  mm-1  of  unit  volume  placed 
at  position  (0,  0,  2)  mm.  The  time-resolved  pro¬ 
files  of  -<57(r,  2,  t)  at  position  r  =  (0,  0,  5)  mm  and 
S7(r,  -2,  t)  at  position  r  =  (0,  2,  0)  mm  from  Eq.  (6) 
and  those  from  diffusion  models  are  shown  in  Fig.  2. 
The  significant  difference  between  our  model  and 
diffusion  models  shows  that  the  nondiffusive  nature 
of  photon  migration  at  early  times  is  important  and 
cannot  be  neglected  when  the  separation  of  any  pair  of 
source,  inhomogeneity,  and  detector  is  small. 

In  conclusion,  a  photon-transport  forward  model  of 
image  reconstruction  in  turbid  media  has  been  derived 
by  use  of  the  Born  approximation  of  th  radiative  trans¬ 
fer  equation.  A  simplified  cumulant  solution  in  an 
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(a) 


t(ps) 

(b) 

Fig.  1.  Photon-distribution  profiles  (a)  J(0)(r,  z,  t)  at  posi¬ 
tion  r  =  (0,  0,5)  mm  and  (b)  7(0)(r,  -z,  t)  at  position  r  = 
(0,  2,  0)  mm. 


(a) 


t(ps) 

(b) 

Fig.  2.  Changes  of  photon-distribution  function  (a)  -8I(r, 
z ,  t)  at  position  r  =  (0,  0,  5)  mm  and  (b)  dl( r,  — z ,  t)  at 
position  r  —  (0,  2,  0)  mm. 


infinite  uniform  medium  serves  as  the  background 
Green’s  function.  This  model  provides  a  correct  pic¬ 
ture  of  nearly  ballistic  motion  of  photons  at  early  times 
and  reduces  to  the  center-moved  diffusion  approxima¬ 
tion  at  long  times.  Extension  to  semi-infinite  and  slab 
geometries  of  this  model  is  being  studied. 
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Parotid  Gland  Tissues  Investigated  by  Picosecond 
Time-Gated  and  Optical  Spectroscopic 
Imaging  Techniques 
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Abstract — Near-infrared  (NIR)  time-resolved  and  spectro¬ 
scopic  transillumination  imaging  techniques  are  used  to  investigate 
normal  tissues  and  Warthin’s  tumor  of  human  parotid  glands. 
The  time-sliced  imaging  arrangement  uses  120-fs,  1-kHz  repeti¬ 
tion-rate,  800-nm  pulses  from  a  Ti :  sapphire  laser  and  amplifier 
system  for  sample  illumination  and  an  ultrafast  gated  intensified 
camera  system  (UGICS)  for  recording  two-dimensional  (2-D) 
images  using  transmitted  light.  Images  recorded  with  earlier 
temporal  slices  (approximately  first  100  ps)  of  transmitted  light 
highlight  the  tumor,  while  those  recorded  with  later  temporal 
slices  (later  than  200  ps)  accentuate  normal  tissues.  The  spectro¬ 
scopic  imaging  arrangement  uses  1210-1300  nm  tunable  output 
of  a  Cr :  forsterite  laser  for  sample  illumination,  a  Fourier  space 
gate  to  discriminate  against  multiple-scattered  light,  and  a  NIR 
area  camera  to  record  2-D  images.  The  tumor  region  in  the 
specimen  appears  brighter  than  the  normal  region  in  spectro¬ 
scopic  images  recorded  with  light  of  different  wavelengths.  A 
wavelength-dependent  variation  in  the  ratio  of  light  transmission 
through  the  tumor  to  that  through  the  normal  parotid  gland  is 
observed.  Differences  in  scattering  and  wavelength-dependent 
absorption  characteristics  of  normal  parotid  gland  and  Warthin’s 
tumor  provide  a  consistent  explanation  of  these  observed  features. 
Histopathological  analysis  of  the  specimen  sheds  light  on  the 
probable  origin  of  the  differences  in  scattering  and  absorption 
characteristics. 

Index  Terms — Near-infrared  imaging,  optical  biomedical 
imaging,  parotid  gland,  photon  propagation  in  highly  scattering 
media,  time-gated  imaging,  transillumination  imaging,  Warthin’s 
tumor. 


I.  Introduction 

OPTICAL  biomedical  imaging  techniques  are  attracting  in¬ 
terest  as  promising  noninvasive  means  for  detection  and 
diagnosis  of  tumors  and  abnormalities  in  human  body  [1]-[1 1]. 
Body  organs  that  are  potentially  amenable  to  optical  imaging  in¬ 
vestigation  include  breast,  brain,  gastro-intestinal  tract,  obstetric 
and  gynecological  tract,  prostate,  skin,  salivary  glands,  arteries, 
and  bones.  Development  of  appropriate  and  effective  optical 
modalities  for  in  vivo  clinical  imaging  of  lesions  in  any  of  these 
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organs  is  an  involved  process.  An  important  step  in  this  devel¬ 
opment  process  is  testing  of  the  efficacy  of  the  method  on  ex 
vivo  tissue  specimens  of  the  target  organ,  which  in  turn  provides 
key  information  about  light  transport  and  optical  spectroscopic 
characteristics  of  the  types  of  tissues  under  investigation.  We 
have  initiated  such  a  study  of  ex  vivo  normal  and  abnormal  tis¬ 
sues  from  different  organs  using  a  variety  of  near-infrared  (NIR) 
optical  imaging  techniques.  In  this  paper,  we  present  results 
of  time-sliced  and  spectroscopic  two-dimensional  (2-D)  NIR 
imaging  experiments  on  excised  normal  tissues  and  Warthin’s 
tumor  of  the  parotid  gland. 

Time-sliced  imaging  [6],  [10],  [12]  is  an  extension  of  the  con¬ 
cept  of  time-resolved  early-light  imaging  [13]— [17]  wherein  a 
sequence  of  2-D  images  of  the  sample  under  investigation  is 
recorded  using  different  temporal  slices  of  light  emergent  from 
the  sample,  in  addition  to  the  early-arriving  part.  Differences  in 
light  transport  properties  of  constituent  tissues  in  a  sample  are 
highlighted  in  the  2-D  images  obtained  with  different  temporal 
slices  of  the  transmitted  light.  We  use  ultrashort  pulses  of  light 
from  a  femtosecond  Tirsapphire  laser  to  illuminate  the  sample 
and  a  time-gated  camera  with  a  variable  gate  position  to  record 
2-D  images  with  the  gated  fraction  of  the  forward  transmitted 
light. 

Spectroscopic  imaging  [18]  uses  light  of  different  colors 
to  exploit  any  spectroscopic  difference  for  enhancing  image 
contrast  and  exploring  diagnostic  potential.  In  a  spectroscopic 
imaging  measurement,  one  uses  light  of  different  wavelengths 
to  record  2-D  transillumination  images  of  the  sample  under 
investigation.  If  a  part  of  the  sample  happens  to  have  a  different 
response  to  light  at  a  particular  wavelength  than  the  rest,  then 
“resonant”  (or  “near-resonant”)  images  recorded  with  light  of 
that  wavelength  (or  a  wavelength  near  the  resonance)  provide 
contrast  between  that  part  and  rest  of  the  sample  [18].  Compar¬ 
isons  with  “nonresonant”  images  that  are  obtained  using  light 
of  wavelengths  away  from  the  resonance  help  accentuate  the 
contrast  even  further. 

The  parotid  glands,  located  in  the  tissue  inferior  and  ante¬ 
rior  to  the  ears,  are  the  largest  of  the  three  main  pairs  of  sali¬ 
vary  glands  [19].  They  secrete  the  serous  (thin,  watery)  com¬ 
ponent  of  the  saliva,  and  the  parotid  duct  delivers  the  secretion 
into  the  mouth.  Diseases  of  the  parotid  glands  include  Warthin’s 
tumor  ( cystadenoma  lymphomatosum)  [20],  and  pleomorphic 
adenoma.  In  this  article,  we  present  results  of  optical  imaging 
experiments  on  an  ex  vivo  tissue  specimen  with  normal  parotid 
gland  tissue  and  Warthin’s  tumor.  Warthin’s  tumor  is  the  second 
most  common  benign  tumor  that  accounts  for  5%  to  10%  of 
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all  parotid  gland  tumors.  The  tumor  commonly  is  a  round  to 
oval  encapsulated  mass  1  cm  to  10  cm  in  diameter  with  a  char¬ 
acteristic  lymphoid  component.  It  is  located  in  the  parotid  or 
preparotid  lymph  nodes  in  99%  of  the  cases  [21],  [22].  The  cut 
surface  of  the  tumor  is  gray  in  color  and  has  cystic  or  cleftlike 
spaces  filled  with  a  gray  or  light  brown  mucinous  secretion.  Tall 
columnar  cells  that  overlie  lymphoid  tissue  line  the  cystic  or 
cleft  like  spaces. 

The  remainder  of  the  article  is  organized  as  follows.  Sec¬ 
tion  II  presents  an  outline  of  the  experimental  methods,  and  a 
brief  description  of  the  samples  used  in  the  experiment.  Sec¬ 
tion  m  presents  optical  imaging  results  along  with  histopatho- 
logical  analysis  of  the  samples.  In  Section  IV,  we  discuss  the 
implications  of  the  experimental  results. 

II.  Experimental  Methods 

A.  Time-Sliced  Imaging 

Experimental  arrangement  for  time-sliced  imaging,  shown 
schematically  in  Fig.  1(a),  made  use  of  approximately  120-fs 
duration,  1-kHz  repetition-rate  pulses  from  a  Ti:  sapphire 
laser  and  amplifier  system  [23]  for  sample  illumination,  and 
an  ultrafast  gated  intensified  camera  system  (UGICS)  for 
recording  2-D  images  using  picosecond-duration  slices  of 
light  transmitted  through  the  sample.  The  laser  output  was 
tunable  over  the  750-850  nm  spectral  range  with  an  average 
power  of  3  W.  We  used  approximately  200  mW  at  800  nm 
for  experiments  reported  here.  A  beam  expander  expanded  the 
beam  and  an  aperture  selected  a  3-cm  diameter  central  part  of 
it  to  illuminate  the  sample. 

The  UGICS  was  a  compact  time-gated  image  intensifier 
unit  that  was  fiber-optically  coupled  to  a  charge-coupled 
device  (CCD)  camera.  It  provided  an  electronic  time  gate 
whose  FWHM  duration  could  be  adjusted  to  a  minimum  of 
approximately  80  ps.  The  temporal  position  of  the  gate  could 
be  electronically  varied  over  a  20-ns  range  with  a  minimum 
step  size  of  25  ps.  The  CCD  camera  had  a  384  x  286  pixels 
sensing  element  with  a  pixel  size  of  23  pm  x  23  pm.  A  24-mm 
focal-length  // 2.8  camera  lens  collected  the  signal  transmitted 
through  the  tissue  sample  and  directed  it  to  the  sensing  element. 
The  transillumination  signal  recorded  by  the  UGICS  at  a 
particular  gate  position  was  a  convolution  of  the  transmitted 
light  pulse  with  the  gate  pulse  centered  on  the  gate  position. 
A  personal  computer  stored  and  displayed  the  resulting  2-D 
images. 

B.  Spectroscopic  Imaging 

The  experimental  arrangement  for  near-infrared  (NIR)  spec¬ 
troscopic  imaging,  displayed  in  Fig.  1(b),  used  the  121Q-1300 
nm  continuous-wave  mode-locked  output  of  a  Cr4+  :  forsterite 
laser  to  illuminate  the  sample.  A  set  of  calibrated  neutral  density 
filters  helped  maintain  the  average  optical  power  of  the  incident 
beam  at  approximately  35  mW  for  all  the  wavelengths  used  in 
the  imaging  experiment.  A  beam  expander  expanded  the  beam 
and  an  aperture  selected  a  3-cm  diameter  central  part  of  it  to  il¬ 
luminate  the  sample. 

A  Fourier  space  gate  [24]  selected  out  a  fraction  of  the  less- 
scattered  image-bearing  photons  from  the  strong  background 
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Fig.  1.  Schematic  diagram  of  the  experimental  arrangement  for:  (a) 
time-sliced  imaging  and  (b)  spectroscopic  imaging.  Key:  A  =  aperture, 
L=  lens,  M  =  mirror,  CCD  =  charge -coupled  device,  UGICS  =  ultrafast  gated 
intensified  camera  system. 


of  the  image-blurring  diffuse  photons.  A  50  mm  focal-length 
camera  lens  placed  on  the  optical  axis  at  a  distance  of  50  mm 
from  the  aperture  in  the  Fourier  gate  collected  and  collimated  the 
low-spatial-frequency  light  filtered  by  the  aperture  and  directed 
it  to  the  128  x  128  pixels  sensing  element  of  an  InGaAs  near- 
infrared  (NIR)  area  camera. 

C.  Sample 

Excised  parotid  gland  tissue  specimen  used  in  the  experi¬ 
ments  reported  in  this  article  was  obtained  from  the  New  York 
Eye  and  Ear  Infirmary  (NYEEI)  under  Internal  Review  Board 
approvals  at  the  City  College  of  New  York  and  NYEEI.  The 
specimen  came  from  the  left  parotid  gland  of  a  37-year-old  male 
patient.  It  was  a  gray/brown  nodule  surrounded  by  normal  sali¬ 
vary  gland  parenchyma.  The  nodule  itself  was  focally  cystic 
with  a  papillary  appearance.  The  clinical  diagnosis  was  a  be¬ 
nign  Warthin’s  tumor  of  the  left  parotid  gland. 

The  specimen  was  received  on  ice.  It  was  held  between  two 
glass  plates  and  slightly  compressed  to  ensure  same  overall 
thickness  for  optical  imaging  measurements.  The  dimensions  of 
this  slightly  compressed  sample  were  20  mm  x  10  mm  x  5  mm. 
Light  transited  through  the  5  mm  path  length  of  the  sample  for 
transillumination  imaging  measurements  reported  here.  After 
NIR  imaging  experiments,  the  sample  was  placed  in  formalin 
and  transferred  to  NYEEI  for  histological  analysis. 
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Fig.  2.  Time-sliced  2-D  images  (left  frames)  of  a  20  mm  x  10  mm  x  5  mm 
parotid  gland  tissue  sample  with  normal  parotid  gland  (N)  and  Warthin’s  tumor 
(T)  for  gate  positions  of  (a)  50  ps,  (b)  100  ps,  and  (c)  250  ps.  Corresponding 
spatial  intensity  profiles  integrated  over  the  same  horizontal  area  highlighted 
by  white  dashed  boxes  in  the  images  are  shown  in  the  right  side  frames.  The 
zero  position  was  taken  to  be  the  time  of  arrival  of  the  light  pulse  through  a 
5-mm-thick  quartz  cell  filled  with  water. 


III.  Results 


Fig.  3.  Spectroscopic  2-D  images  (left  frames)  of  the  parotid  gland  tissue 
sample  with  normal  parotid  gland  (N)  and  Warthin’s  tumor  (T)  for  light  of 
wavelengths  (a)  1225  nm,  (b)  1250  nm,  and  (c)  1285  nm.  Corresponding  spatial 
intensity  profiles  integrated  over  the  same  horizontal  area  highlighted  by  white 
dashed  boxes  in  the  images  are  shown  in  the  right  frames.  No  time  gating  was 
used,  (d)  Spatial  profiles  of  the  ratio  of  intensities  R$p(A),  as  described  in  the 
text,  for  the  reference  wavelength  of  1285  nm,  and  probing  wavelengths  of  1225 
and  1250  nm. 


A.  Time-Sliced  Imaging 

Time-sliced  2-D  transillumination  images  of  the  parotid 
gland  sample  for  gate  positions  of  50  ps,  100  ps,  and  250  ps 
are  displayed  in  the  left  frames  of  Fig.  2(a)-(c),  respectively. 
The  zero  position  was  taken  to  be  the  time  of  arrival  of  the  light 
pulse  through  a  5-mm-thick  quartz  cell  filled  with  water.  The 
corresponding  frames  to  the  right  present  the  spatial  intensity 
profiles  of  the  respective  images  integrated  over  the  same  hori¬ 
zontal  area  in  all  the  figures.  The  salient  feature  of  the  images 
is  the  differences  in  time-dependent  brightness  of  the  normal 
tissue  and  Warthin’s  tumor  in  the  sample.  In  the  early  50-ps 
image  of  Fig.  2(a)  only  the  tumor  is  visible.  The  corresponding 
horizontal  spatial  profile  exhibits  peak  in  the  tumor  region  and 
hits  the  baseline  in  the  normal  tissue  region.  At  this  early  time, 
markedly  more  light  is  transmitted  through  the  tumor  than  the 
normal  parotid  gland  tissue.  With  time  the  image  of  normal 
parotid  gland  tissue  gains  in  brightness  and  the  corresponding 
region  of  the  spatial  profile  rises  above  the  baseline,  as  the 
typical  sampling  of  Fig.  2(b)  for  the  delay  position  of  100  ps 
shows.  At  later  times  (250  ps)  the  image  of  the  normal  gland 
becomes  brighter  than  the  tumor,  the  corresponding  spatial 


profile  peaks  in  the  normal  gland  position  indicating  higher 
light  transmission  through  the  normal  gland  than  the  tumor 
[Fig.  2(c)]. 

These  results  indicate  that  light  transits  faster  through 
Warthin’s  tumor  than  normal  parotid  gland  tissue.  Lower 
scattering  or/and  higher  absorption  of  light  by  the  tumor  may 
account  for  the  observed  temporal  behavior.  Since  there  is  no 
known  absorption  of  800-nm  light  by  parotid  gland  tissues,  we 
attribute  these  time-dependent  differences  in  the  relative  light 
transmission  through  Warthin’s  tumor  and  normal  parotid  gland 
tissues  to  the  differences  in  the  light  scattering  characteristics 
of  these  components  of  tissues. 

B.  Spectroscopic  Imaging 

Spectroscopic  images  of  the  specimen  recorded  using  light 
of  wavelengths  1225, 1250,  and  1285  nm  appear  in  the  left  side 
frames  of  Fig.  3(aMc),  respectively.  Corresponding  spatial  pro¬ 
files  integrated  over  the  same  horizontal  area  in  all  the  images 
appear  in  the  right  side  frames.  The  salient  feature  of  the  images 
is  the  higher  brightness  of  the  Warthin’s  tumor  region  compared 
to  the  normal  parotid  gland  region  for  all  three  wavelengths. 
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Since  no  time  gating  is  used  in  spectroscopic  imaging,  the  re¬ 
sult  is  indicative  of  higher  overall  light  transmission  through  the 
tumor  than  that  through  the  normal  parotid  gland. 

Another  subtle  feature  is  the  wavelength  dependence  of 
relative  light  transmission  through  the  tumor  and  the  normal 
tissue.  We  examined  the  variation  in  relative  brightness  with 
wavelength  by  monitoring  the  tumor-to-normal  intensity 

ratio,  JZtm(A)  7tumor(A)//n0rmal(A),  where  /tumor(A) 

/normal (A),  are  the  averaged  intensities  through  the  Warthin’s 
tumor  and  normal  parotid  gland  regions,  respectively,  of  the 
spatial  profile  at  wavelength  A.  The  values  of  J?TN(A)  for  1225, 
1250,  and  1285  nm  are  2.9,  2.0,  and  1 .7,  respectively. 

Even  more  instructive  is  the  ratio  of  intensities  i?SP(A)  = 
7(A)//(Ar),  where  /(A)  and  /( Ar),  are  the  values  of  inten¬ 
sities  at  probing  wavelength  A  and  selected  reference  wave¬ 
length  Ar,  respectively,  measured  at  the  same  location  in  the 
respective  images.  The  horizontal  spatial  profiles  of /?sp(A)  for 
Ar  =  1285  nm  and  A  =  1225  and  1250  nm  averaged  over 
the  width  of  the  dashed  lines  in  the  images  of  Fig.  3(a)— (c)  are 
shown  in  Fig.  3(d).  These  profiles  were  obtained  by  taking  the 
ratio  of  spatial  profile  at  wavelength  A  to  that  at  1285  nm.  The 
salient  feature  of  the  i?sp(A)  profiles  in  Fig.  3(d)  is  that  the 
tumor  region  tends  to  have  a  higher  value  of  the  ratio  (>  1 )  than 
the  normal  region  (<!)• 

The  transition  region  between  the  Warthin’s  tumor  and 
normal  parotid  gland  appears  brighter  than  the  rest  of  the  tumor 
in  the  spectroscopic  images  and  the  corresponding  spatial 
intensity  profiles  tend  to  peak  around  pixel  50.  The  peak  is 
more  prominent  in  the  spatial  profile  of  Fig.  3(c)  than  in  other 
profiles.  A  histological  examination  of  the  specimen  revealed 
that  although  this  region  was  composed  of  the  same  lymphoid 
and  epithelial  elements  as  the  rest  of  the  tumor,  the  tissue  in  this 
area  was  less  dense  and  had  more  lymphoid  tissue  than  epithe¬ 
lial  tissue  compared  to  the  rest  of  the  tumor.  These  differences 
lead  to  higher  light  transmission  through  the  transition  region. 
This  result  indicates  the  sensitivity  of  spectroscopic  imaging  to 
variations  in  tissue  density  and  morphology. 

IV.  Discussion 

The  results  of  both  the  time-sliced  imaging  and  spectroscopic 
imaging  experiments  show  the  ability  of  the  two  techniques  to 
select  between  the  Warthin’s  tumor  and  normal  parotid  gland 
tissues.  The  contrast  between  the  Warthin’s  tumor  and  the 
normal  parotid  gland  is  most  dramatic  in  the  images  recorded 
with  earlier  and  later  time  slices.  We  observed  similar  contrast 
between  normal  and  cancerous  regions  in  the  time-sliced 
images  of  human  breast  tissue  specimens  as  well  [6],  [25], 

Spectroscopic  images  exhibit  three  main  features:  (a)  overall 
higher  light  transmission  through  the  tumor;  (b)  a  wave¬ 
length-dependent  variation  in  tumor-to-normal  intensity  ratio, 
#tn(A);  and  (c)  a  variation  in  the  spatial  profile  of  the  ratio, 
iisp(A)  between  the  tumor  and  normal  regions  of  a  specimen. 
We  attribute  the  higher  light  transmission  through  the  Warthin’s 
tumor  compared  to  that  through  the  normal  parotid  gland  at 
all  three  wavelengths  to  lower  light  scattering  by  the  tumor. 
This  attribution  is  consistent  with  the  results  of  the  time-sliced 
imaging  experiment  carried  out  with  800-nm  femtosecond 


(b) 


Fig.  4.  Histological  micrograph  of  a  representative  section  from:  (a)  the 
normal  parotid  gland  used  in  the  optical  imaging  experiment  showing  fat  cells 
(FC),  glandular  alveoli  (GA),  and  ducts  (Duct)  and  (b)  the  Warthin’s  tumor 
showing  double  layers  (columnar  and  cuboidal)  of  ductal  epithelium  (DE), 
clefts  of  the  duct  (CD),  papillary  projections  (PP)  of  the  ductal  epithelium, 
cystic  spaces  (CS),  and  areas  of  lymphoid  tissue  (LT). 

light  pulses.  The  wavelength  dependent  variation  in  f?Tx(A) 
and  iZsp(A)  is  expected  to  be  associated  with  difference  in 
absorption  by  the  normal  parotid  gland  and  tumor. 

In  order  to  obtain  a  better  appreciation  of  the  origin  of  the  dif¬ 
ferences  in  scattering  and  absorption  characteristics  of  normal 
parotid  gland  and  Warthin’s  tumor,  we  turned  to  the  histopatho- 
logical  analysis  of  the  sample.  Histopathology  provides  detailed 
information  about  the  structural  differences  between  the  normal 
and  abnormal  tissues  and  sets  the  standard  against  which  the 
performance  of  any  new  technique  should  be  evaluated. 

Fig.  4(a)  presents  a  histological  micrograph  of  a  repre¬ 
sentative  section  of  the  normal  parotid  gland  tissue  used  in 
the  optical  imaging  experiment.  The  key  features  are  the  fat 
cells  (FC),  glandular  alveoli  (GA),  and  ducts  (Duct).  Fig.  4(b) 
is  a  micrograph  of  a  similar  representative  section  from  the 
Warthin’s  tumor  in  the  specimen  under  same  magnification. 
The  features  in  the  histology  of  the  Warthin’s  tumor  are  double 
layers  (columnar  and  cuboidal)  of  ductal  epithelium  (DE), 
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clefts  of  the  duct  (CD),  papillary  projections  (PP)  of  the  ductal 
epithelium,  cystic  spaces  (CS),  and  areas  of  lymphoid  tissue 
(LT).  The  contrast  with  the  histology  of  the  normal  parotid 
gland  is  quite  distinct. 

Light  scattering  depends  on  the  size,  shape,  and  relative  re¬ 
fractive  indices  of  the  scattering  entities.  The  aforementioned 
structural  differences  between  the  tumor  and  normal  parotid 
gland  are  likely  to  lead  to  different  light  scattering  and  spectro¬ 
scopic  characteristics.  A  quantitative  evaluation  of  the  role  of 
different  scatterers  is  highly  involved  that  will  require  indepen¬ 
dent  measurement  of  the  optical  and  scattering  characteristics 
of  each  of  the  different  scattering  entities,  and  modeling  of  their 
distribution  in  the  tissue  ultrastructure. 

A  particularly  significant  difference  observed  in  the  micro¬ 
graphs  of  the  two  representative  sections  that  has  important  con¬ 
sequences  for  NIR  imaging  experiments  is  the  paucity  of  fat 
(adipose)  cells  in  the  Warthin’s  tumor,  and  their  abundance  in 
the  normal  parotid  gland.  Two  characteristics  of  adipose  tissues 
may  shed  light  on  the  observed  results  of  NIR  imaging  experi¬ 
ments.  First,  adipose  tissues  scatter  light  more  effectively  than 
other  types  of  tissues,  such  as,  fibrous,  glandular,  and  ductal 
[6],  [25],  [26].  Second,  adipose  tissue  has  a  broad  optical  ab¬ 
sorption  band  centered  on  1203  nm  [26].  The  histological  anal¬ 
ysis  shows  that  the  normal  parotid  gland  tissue  has  a  much 
higher  abundance  of  adipose  cells  than  that  in  Warthin’s  tumor. 
Light  transiting  through  the  normal  tissue  will  be  more  scat¬ 
tered  by  the  adipose  cells,  traverse  longer  distance  within  the 
sample,  and  arrive  later  in  time  than  that  transiting  through  the 
Warthin’s  tumor.  It  thus  provides  an  explanation  of  the  results  of 
time-sliced  imaging  experiments  that  show  early-light  images 
highlighting  the  tumor,  while  the  late-light  images  accentuate 
the  normal  parotid  gland. 

The  results  of  spectroscopic  imaging  arrangement  may  be  ex¬ 
plained  as  well.  The  normal  part  of  the  sample  appears  darker 
in  the  spectroscopic  images  Fig.  3(a)-(c) as  the  more  abundant 
adipose  cells  there  contribute  to  higher  scattering  (more  loss  of 
light)  than  the  tissues  in  the  tumor.  Adipose  cells  in  the  normal 
parotid  gland  are  one  of  the  major  contributors  to  the  scattering, 
but  may  not  be  the  sole  reason  for  higher  scattering  by  normal 
tissue. 

The  light  wavelength  of  1225  nm  is  closer  to  the  adipose  ab¬ 
sorption  peak  at  1203  nm,  1250  nm  is  in  the  wing  of  the  ab¬ 
sorption  band,  while  1285  nm  is  away  from  the  absorption  reso¬ 
nance.  Consequently,  among  the  three  wavelengths  presented  in 
Fig.  3(a)-(c),  1225-nm  light  is  absorbed  the  most  by  the  adipose 
cells  in  the  normal  part  of  the  specimen,  1285-nm  light  the  least, 
with  the  1250-nm  light  in  between.  This  wavelength-dependent 
absorption  contributes  to  the  observed  decrease  in  the  ratio  be¬ 
tween  1225  nm  to  1285  nm.  The  tumor-to-normal  ratio  jRtn(A) 
thus  acts  as  a  parameter  to  differentiate  between  normal  tissue 
and  Warthin’s  tumor  of  the  parotid  gland.  The  use  of  this  ratio 
for  in  vivo  application  requires  knowledge  of  the  suspect  region 
and  a  normal  region,  and  evaluation  of  the  ratio  between  trans¬ 
mitted  intensities  through  the  abnormal  to  normal  region. 

The  ratio  Rsp{ty  is  a  more  promising  parameter  since  it  in¬ 
volves  measurements  of  two  intensities,  one  at  a  reference  wave¬ 
length  and  the  other  at  a  wavelength  near-resonant  with  adipose 
absorption,  at  the  same  location  and  evaluation  of  the  ratio.  In 


obtaining  the  spatial  profiles  of  R$p(\)  displayed  in  Fig.  3(d), 
we  used  1285  nm  that  is  far-removed  from  the  adipose  reso¬ 
nance  peak  at  1203  nm  as  the  reference  wavelength,  and  1225 
and  1250  nm  that  lie  within  the  adipose  optical  absorption  band 
as  probing  wavelengths.  We  obtain  a  higher  contrast  between 
the  tumor  and  normal  regions  in  the  RspW  profile  that  uses 
intensities  at  1225  nm,  the  wavelength  closer  to  the  resonance 
peak.  However,  in  both  the  i?sp(A)  profiles  using  1225  nm  and 
1250  nm,  the  ratio  values  tend  to  be  more  than  one  for  the  tumor 
region,  and  less  than  1  for  the  normal  region.  Assuming  that  this 
trend  holds,  a  measurement  of  1?sp(A)  value  at  a  location  in  the 
specimen  may  identify  it  as  either  tumor  (R$ p  >  1)  or  normal 
{Rsp  <  1)»  and  no  comparison  with  another  region  is  needed. 
More  measurements  on  a  larger  number  of  samples  with  dif¬ 
ferent  stages  of  tumor  progression  will  be  needed  to  build  up 
statistics  for  establishing  the  values  of  parameters  R$p  and  Rt  n 
that  will  be  indicative  of  the  status  of  tissue  (tumor  or  normal) 
with  requisite  specificity. 

In  summary,  the  results  of  these  experiments  further  demon¬ 
strate  the  potential  of  NIR  time-sliced  and  spectroscopic 
imaging  techniques  for  detection  and  diagnosis  of  tumors  in 
optically  accessible  suspect  body  organs. 
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The  recent  controversy  over  whether  the  photon  diffusion  coefficient  depends  on  absorption  is  addressed 
by  use  of  the  analytical  solution  of  the  photon  transport  equation  in  an  infinite  homogeneous  scattering 
medium.  The  diffusion  coefficient  is  found  to  be  independent  of  absorption  but  temporally  dependent.  Af¬ 
ter  a  long  period  of  time,  the  photon  diffusion  coefficient  approaches  D  =  1/3 /V,  which  supports  a  claim 
made  by  Furutsu  and  Yamada  [Phys.  Rev.  E  50,  3634  (1994)].  At  early  times,  the  diffusion  coefficient  is 
smaller  than  D  *  1/3/i/,  but  this  reduction  cannot  be  expressed  as  D  =  l/3(/x/  +  pa),  since  the  time- 
dependent  diffusion  coefficient  is  found  to  be  unrelated  to  absorption.  ©  2002  Optical  Society  of  America 
OCIS  codes:  290.1990,  030.5620,  170.5280,  290.7050. 


The  diffusion  approximation  of  the  radiative  trans¬ 
fer  equation  has  been  used  extensively  in  the  study 
of  photon  migration  and  applications  in  turbid  media. 
A  controversy  exists  about  the  form  of  diffusion  coef¬ 
ficient  D  in  media  in  coexistence  with  scattering  and 
absorption.  In  traditional  derivations  of  the  diffusion 
equation  from  the  transport  equation,  the  diffusion  co¬ 
efficient  is  obtained  as1-4 

D-1/3(ai.'  +  m«).  (1) 

where  fial  is  the  reduced  scattering  coefficient  and  jia 
is  the  absorption  coefficient.  Furutsu  and  Yamada6 
and  Furutsu6  pointed  out  that  this  form  of  D  does 
not  consist  of  the  transformation  property  of  the 
transport  equation  in  an  infinite  homogeneous  tur¬ 
bid  medium:  the  photon  distribution  function  at 
position  r,  in  direction  s  and  time  t ,  /( r,  s,  £)  — 
exp(-cjjLat)Io(r ,  s,  t),  where  /o(r,  s,  *)  is  the  solution 
for  the  same  medium  but  without  absorption.  They5,6 
suggested  another  derivation  that  leads  to 

D  -  lt/ 3  -  1/3/ttV,  (2) 

which  is  independent  of  absorption,  where  lt  is  the 
transport  mean-free  path.  This  claim  has  been  sup¬ 
ported  by  Bassani  et  al. ,7  Durduran  et  a/.,8  and  Nakai 
et  al9  based  on  experiments  and  Monte  Carlo  simula¬ 
tions.  In  contrast,  in  later  papers  Aronson  and  Corn- 
gold,10  Rinzema  et  al.,11  and  Durian12,  based  on  their 
experiments  and  numerical  simulations,  asserted  that 
D  should  depend  on  absorption. 

Most  recently,  we  developed  an  analytical  solution 
of  the  transport  equation  in  an  infinite  homogeneous 
medium.13,14  We  derived14  an  expression  for  the 
spatial  cumulants  of  the  photon  distribution  function 
7(r,  s,  t)  at  any  direction  s  and  time  t ,  exact  up  to  an 
arbitrarily  high  cumulant  order,  which  can  be  used  for 
accurate  computation  of  the  photon  distribution  func¬ 
tion.  Up  to  the  second  cumulant  order,  we  obtained 
an  approximate  Gaussian  spatial  distribution  that 
has  the  exact  central  position  and  the  exact  half-width 
of  the  distribution.13  For  photon  density  N( r,  t)y 
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at  t  — ►  oo  our  result  approaches  the  center-moved 
diffusion  approximation  (CMDA),  with  D  —  1/3 pig 
independent  of  absorption.  At  a  finite  time,  the 
diffusion  coefficient  D  increased  from  zero  at  t  =  0 
to  the  above  value  at  t  — ►  ®.  Typically,  in  the  case  of 
the  anisotropic  factor  g  ~  0.9  it  takes  approximately 
10Z*  for  the  diffusion  mode  to  be  valid,  as  shown 
experimentally  by  Yoo  et  al.15  This  D(t)  as  a  function 
of  time  is  determined  only  by  scattering  parameters 
unrelated  to  absorption.  Hence,  it  is  physically  un¬ 
reasonable  to  use  Eq.  (1)  to  reduce  the  value  of  D  at  a 
finite  time  to  fit  experimental  or  simulated  data. 

The  time-dependent  photon  density  N( r,  t)  obtained 
from  the  CMDA  in  an  infinite  homogeneous  medium 
for  a  collimated  pulse  source  located  at  r  —  0  with  in¬ 
cident  direction  along  z  is  given  by16 

"<'•‘>-^+^--4 <3) 

where  D  is  the  diffusion  constant.  This  solution 
differs  from  the  standard  diffusion  solution  in  that  the 
center  of  the  distribution  is  moved  by  lit  from  the 
source  position  along  the  incident  direction.  On 
the  other  hand,  our  cumulant  approximation  (CUMA) 
for  photon  density,  exact  up  to  the  second-order 
cumulant  from  the  same  source,  is  obtained  as13,14 

1 _ 1_  f  iz-Rz)2 1 

N{T'  l)  (4n D„ct)W  4irDxxct  eXP|_  4 D„ct  J 

x  exp[ "  ]exp(-^c*)  >  (4) 

with  the  moving  center  located  at 

Rz  =  c[l  -  exp(-gi*)]/gi .  (5) 

The  corresponding  diffusion  coefficients  are  given  by 
©  2002  Optical  Society  of  America 
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£«  = 


t_  _  Sgi  -  g2 
gl  gl2(gl  ~  gi) 


[1  -  exp(-gif)] 


giigi  ~  gi) 


[1  -  exp(-g2*)] 


2*i2 


[1  -  exp(-gi<)]2 


(6) 


Dxx  Dyy  3 1 


i+e,C-8i)cl~exp(^it>1 


1 

g*(gl  ~  g2) 


[1  -  exp(-g2*)] 


(7) 


Here  gi  =  /i*c[  1  -  ai/(2l  +  1)],  where  the  single¬ 
scattering  phase  function  is  expanded  in  Legendre 
polynomials  by 


P(cos  0)  =  [l/(47r)]  ^aiPi( cos  0) . 
i 

Two  special  values  of  gi  are  go  =  0  and  gi  =  c/lt- 

The  original  meaning  of  diffusion  is  a  description  of 
a  random  process.  As  early  as  when  Brownian  motion 
was  first  studied,  researchers  have  known  that  under 
random  forces  from  the  surrounding  medium  particles 
will  take  a  diffusion  process,  namely,  spread  from  the 
center  outward.  The  diffusion  coefficient  is  a  parame¬ 
ter  that  characterizes  the  rate  of  spread,  which  can 
be  time  dependent.  The  standard  diffusion  equation 
with  a  diffusion  constant  is  the  simplest  form  with 
which  to  describe  these  phenomena.  The  photon 
propagation  in  a  turbid  medium  is  one  example  of  the 
random  processes,  which  is  more  complicated  than  that 
described  by  the  standard  diffusion  equation,  mainly 
because  the  photons  are  injected  with  velocity  along 
a  direction.  This  leads  to  photon  propagation  from 
ballistic  to  snakelike  and  then  to  the  diffusive  mode. 
The  diffusion  coefficient  in  this  process  should  be  time 
dependent. 

The  fundamental  time-independent  parameters 
of  the  medium  are  the  scattering  coefficients,  the 
absorption  coefficients,  and  the  phase  function  in  the 
radiative  transfer  equation,  because  they  have  definite 
meaning  from  a  microscopic  viewpoint.  When  regard¬ 
ing  the  diffusion  constant  as  a  physical  parameter,  we 
should  note  that  this  concept  is  not  original  but  is  a 
secondary  quantity  that  is  derived  by  use  of  an 
approximate  method.  Of  course,  these  two  meanings 
of  a  diffusion  coefficient  are  not  incompatible.  When 
a  standard  diffusion  equation  is  approximately  valid, 
the  diffusion  constant  is  correct  according  to  both 
meanings. 

Although  the  Gaussian  distribution  in  Eq.  (4)  is  an 
approximation  because  it  cuts  into  the  second-order  cu- 
mulant,  Eqs.  (6)  and  (7)  provide  an  exact  description  of 
the  half-width  of  the  real  distribution. 


Figure  1  shows  the  diffusion  coefficients  from 
CUM  A,  Dzz  and  Dxx  [Eqs.  (6)  and  (7)],  as  a  function  of 
time,  where  gi  are  calculated  by  Mie  theory17  assum¬ 
ing  (for  this  figure)  that  the  homogeneous  scattering 
medium  consists  of  water  droplets  with  p/A  =  1 
uniformly  distributed  in  air,  with  p  the  radius  of 
the  droplet,  A  the  wavelength  of  light,  and  index  of 
refraction  m  =  1.33.  The  anisotropic  factor  for  this 
case  is  ai/3  =  0.8436. 

At  t  —>  oo,  the  photon  density  in  Eq.  (4)  approaches 
Eq.  (3),  where  the  diffusion  coefficients  Dzz  and  Dxx 
approach  D  —  1/3 /**',  not  Eq.  (1).  This  means  that 
Eq.  (1)  for  a  diffusion  coefficient  that  depends  on 
absorption  leads  to  an  incorrect  result  even  for  an 
infinite  time  limit.  On  the  other  hand,  the  central- 
limit  theorem  claims  that  the  obtained  Gaussian 
distribution  would  be  accurate  after  a  large  number  of 
collisions.  Hence,  Eq.  (3)  with  D  —  l/S/ia'  provides 
an  accurate  solution  of  photon  transport  in  an  infinite 
homogeneous  medium  at  an  infinite  time  limit,  no  mat¬ 
ter  whether  absorption  is  strong  or  weak.  At  finite 
times,  the  diffusion  coefficient  is  always  smaller  than 
its  value  at  infinite  time.  In  fact,  at  early  time  t  ->  0, 
the  center  moves  as  ctz  and  the  diffusion  coefficient 
approaches  zero.  This  result  presents  a  clear  picture 
of  nearly  ballistic  motion  at  t  — * ►  0.  With  an  increase 
in  time,  the  motion  at  the  center  slows  down,  and 
the  diffusion  coefficient  increases  from  zero.  This 
stage  of  photon  migration  is  in  the  snakelike  mode. 
For  a  large  period  of  time,  Eq.  (4)  approaches  the 
center-moved  diffusion  solution.  As  shown  in  Eqs.  (6) 
and  (7),  the  time-dependent  diffusion  coefficient  is 
determined  only  by  the  scattering  parameters  gi  and 
is  unrelated  to  the  absorption  coefficient  pa-  Hence, 
to  fit  data  of  numerical  simulation  or  experiments  at 
finite  time,  one  should  not  use  Eq.  (1)  to  reduce  the 
diffusion  coefficient. 

The  Aronson  and  Comgold  paper10  emphasizes 
that  the  time-independent  diffusion  equation  rather 
than  the  time-dependent  equation  should  be  examined 
to  determine  which  form  of  diffusion  coefficient  is 
correct.  The  time-independent  solution  can  be  ob¬ 
tained  by  integration  of  the  time-dependent  solution 
over  time.  Integration  of  Eq.  (3)  over  time  yields 
the  time-independent  diffusive  solution  in  an  infinite 
homogeneous  medium: 


Fig.  1.  Diffusion  coefficients  Dzz  and  Dxx  from  Eqs.  (6) 
and  (7)  as  a  function  of  time  t. 
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Fig.  2.  Steady-state  photon  density  as  a  function  of  dis¬ 
tance  from  the  source  (along  the  incident  direction)  for  dif¬ 
ferent  absorption  coefficients  pa  with  unit  l/lt  obtained 
with  the  CMDA  and  the  CUM  A.  The  unit  of  length  is  lt; 
the  unit  of  time  is  lt/c. 


N{r)  ~  4nDc\l  -  l,z\  eXpHr  ~  WP)1/2]-  (8) 

A  more  accurate  steady-state  distribution  N( r)  in  an 
infinite  homogeneous  medium  can  be  obtained  by  inte¬ 
gration  of  Eq.  (4)  over  time  from  t  =  0  to  t  =  Fig¬ 
ure  2  shows  N( r)  as  a  function  of  r  with  the  detector 
set  at  (0,  0,  r)  for  different  absorption  coefficients  pa- 
The  dashed  curves  were  obtained  from  Eq.  (8)  with 
D  =  it! 3,  and  the  solid  curves  were  obtained  by  in¬ 
tegration  of  our  time-dependent  cumulant  solution  in 
Eq.  (4)  over  time,  with  the  corresponding  diffusion  coef¬ 
ficients  given  by  Eqs.  (6)  and  (7).  When  the  detector  is 
near  the  source,  the  photon  density  from  CUMA  is  dis¬ 
tinctly  smaller  than  that  from  Eq.  (8).  This  happens 
because  early  time  photons  play  an  important  role,  and 
the  diffusion  approximation  fails  when  it  is  near  the 
source.  With  an  increase  in  distance  between  source 
and  detector,  the  results  from  the  CMDA  and  CUMA 
approaches  are  in  agreement  for  different  absorption 
coefficients.  These  results  confirm  that  the  diffusion 


coefficient  should  be  D  =  lt/S  at  the  diffusion  limit, 
independently  of  absorption. 

In  conclusion,  our  analytical  solution  of  the  time- 
dependent  photon  density  in  an  infinite  homogeneous 
medium  supports  the  claim  that  the  diffusion  coeffi¬ 
cients  is  independent  of  absorption  but  is  temporally 
dependent. 
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A  photon  transport  model  for  light  migration  in  turbid  media  based  on  a  cumulant  approximation  to  radiative 
transfer  is  presented  for  image  reconstruction  inside  an  infinite  medium  or  a  bounded  medium  with  a  planar 
geometry.  This  model  treats  weak  inhomogeneities  through  a  Bom  approximation  of  the  Boltzmann  radiative 
transfer  equation  and  uses  the  second-order  cumulant  solution  of  photon  density  to  the  Boltzmann  equation  as 
the  Green's  function  for  the  uniform  background.  It  provides  the  correct  behavior  of  photon  migration  at  early 
times  and  reduces  at  long  times  to  the  center-moved  diffusion  approximation.  At  early  times,  it  agrees  much 
better  with  the  result  from  the  Monte  Carlo  simulation  than  the  diffusion  approximation.  Both  approximations 
agree  well  with  the  Monte  Carlo  simulation  at  later  times.  The  weight  function  for  image  reconstruction  under 
this  proposed  model  is  shown  to  have  a  strong  dependence  at  both  early  and  later  times  on  absorption  and/or 
scattering  inhomogeneities  located  in  the  propagation  direction  of  and  close  to  the  source,  or  in  the  field  of 
view  of  and  close  to  the  detector.  This  effect  originates  from  the  initial  ballistic  motion  of  incident  photons, 
which  is  substantially  underestimated  by  the  diffusion  approximation. 

DOI:  10. 1 103/PhysRevE. 65. 066609  PACS  number(s):  42.25.Dd,  05.60.Cd,  42.30. Wb,  02.70.Uu 

I.  INTRODUCTION  port  mean  free  path  /f~lmm  for  human  breast  tissue.  The 

diffusion  approximation  is  thus  incorrect  for  a  substantial 
Photon  migration  in  turbid  media  is  a  random  walk  in  scattering  range.  In  optical  tomography  [8-13]  where  the 

which  rays  or  photons  traverse  a  medium  of  scatterers  and  distribution  of  inhomogeneities  inside  a  highly  scattering 

absorbers,  and  undergo  multiple  scattering  and  absorption  medium  is  reconstructed  from  measurements  of  the  transmit- 

events  before  escaping.  A  natural  framework  to  deal  with  this  ted  light  surrounding  the  medium,  the  diffusion  approxima- 

type  of  problem  is  provided  by  the  theory  of  radiative  trans-  tion  yields  a  much  underestimated  weight  function  when  any 

fer  in  Chandrasekhar’s  classic  text  [1].  The  linear  Boltzmann  separation  between  the  source,  the  inhomogeneity,  and  the 

equation  governs  the  radiation  field  in  a  medium  that  ab-  detector  is  small.  This  error  may  distort  the  signal  from  the 

sorbs,  emits,  and  scatters  radiation  [2].  Because  the  Boltz-  inhomogeneity  inside  the  medium  because  the  weight  func- 

mann  equation  is  a  nonseparable,  integro-differential  equa-  tion  near  surface  is  usually  much  larger  than  that  inside, 

tion  of  first  order  for  which  an  exact  closed-form  solution  is  Recently,  an  analytical  solution  to  the  Boltzmann  equation 

not  known  except  for  a  few  special  cases,  various  approxi-  was  derived  by  the  authors  in  an  infinite  uniform  medium 

mations  have  been  devised  [1,3,4].  The  most  common  ap-  using  a  cumulant  expansion  [14,15].  An  exact  but  formal 

proximation  is  the  diffusion  approximation,  which  corre-  solution  to  the  Boltzmann  equation  yields  the  photon  distri- 

sponds  to  the  lowest-order  truncation  in  the  spherical  bution  function  J(r,s,0  at  position  r,  direction  s,  and  time  f, 

harmonic  expansion  of  the  photon  distribution  function.  It 
follows  from  the  Boltzmann  equation  under  the  assumption 
that  the  photon  distribution  is  almost  isotropic  after  a  suffi¬ 
cient  large  number  of  scattering  events,  and  thus  provides  an 
asymptotic  approximation  applicable  to  later  times  [5].  The  for  a  source  £(r— r0)£(s—  s0)£(0,  where  ()  means  an  en- 

diffusion  approximation  is  invalid  when  the  incident  photon  semble  average  in  photon  direction  space.  Equation  (1)  is 

still  retains  its  directionality  preference.  Moreover,  approxi-  evaluated  in  Fourier  space  with  the  use  of  the  well-known 

mations  using  higher-order  truncation  in  the  spherical  har-  cumulant  expansion  theorem  [16].  An  algebraic  closed  form 

monic  expansion  of  the  photon  distribution  function  are  still  of  expression  is  obtained  for  an  arbitrary  «th  order  cumulant. 

inefficient  in  describing  the  ballistic  movement  of  photons  at  This  expansion  is  inherently  different  from  the  spherical  har- 

early  times  [6].  Yoo  et  al.  [7]  reported  that  the  diffusion  ap-  monies  expansion  of  the  photon  distribution.  The  first-order 

proximation  fails  for  small  and  intermediate  scattering  cumulant  calculation  determines  the  exact  center  position  of 

ranges.  The  range  of  failure  is  proportional  to  the  transport  the  photon  distribution;  the  second-order  cumulant  calcula- 

mean  free  path  lt~lsf(l  ~g)  where  ls  is  the  scattering  mean  tion  determines  the  exact  half  width  of  the  photon  distribu- 

free  path  and  g  is  the  scattering  anisotropy  (the  average  co-  tion  in  addition;  higher-order  cumulant  calculations  provide 

sine  of  the  scattering  angle).  For  one  important  class  of  ap-  progressively  more  details  of  the  shape  of  the  photon  distri- 

plications  of  photon  migration  in  a  turbid  medium — the  bution  but  do  not  modify  the  cumulants  of  lower  order.  This 

medical  applications,  the  medium  has  a  strongly  peaked  is  a  major  advantage  of  the  cumulant  expansion.  The  photon 

phase  function  in  the  forward  direction  and  a  typical  trans-  distribution  approaches  a  Gaussian  distribution  as  the  num¬ 

ber  of  scattering  events  increases  according  to  the  central 
limit  theorem  [16].  So  it  is  not  surprising  that  the  second- 
*Electronic  address:  minxu@sci.ccny.cuny.edu  order  cumulant  solution  with  a  correct  center  position  and 


7(r,s,0=  S\  r 


r— c  I  s (t')dt' 
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half-width  has  already  provided  a  clear  picture  of  the  time 
evolution  of  photon  migration  from  the  initial  ballistic  to  the 
final  diffusive  regime — that  photons  migrate  with  a  center 
that  advances  in  time,  and  with  an  ellipsoidal  contour  that 
grows  and  changes  shape  [14]. 

The  cumulant  solution  depends  explicitly  on  the  phase 
function  of  the  medium  and  involves  a  complicated  numeri¬ 
cal  integration  over  angular  parameters  to  build  a  forward 
model.  It  is  inconvenient  for  direct  use  in  image  reconstruc¬ 
tion.  An  approximate  form  of  the  second-order  cumulant  so¬ 
lution  relating  the  scattered  wave  field  directly  to  the  weak 
inhomogeneities  in  an  infinite  space  was  later  proposed  by 
the  authors  [17],  which  retains  the  main  features  of  photon 
propagation  at  both  early  and  later  times  and  reduces  to  the 
conventional  diffusion  approximation  at  later  times. 

In  this  paper,  we  will  first  extend  the  second-order  cumu¬ 
lant  solution  to  planar  geometries  (semi-infinite  and  slab  me¬ 
dia)  after  a  brief  recount  of  the  main  results  of  the  cumulant 
solution  to  the  Boltzmann  equation  in  an  infinite  space.  The 
result  of  Monte  Carlo  simulations  is  then  presented  for  both 
infinite  and  semi-infinite  media  to  verify  the  behavior  of  the 
second-order  cumulant  solution  at  both  early  and  later  times. 
The  weight  function  for  image  reconstruction  of  weak  inho¬ 
mogeneities  is  calculated  with  use  of  the  simplified  cumulant 
and  diffusion  approximations  for  semi-infinite  and  slab  me¬ 
dia.  The  results  from  the  two  approximations  are  compared. 
The  advantage  of  this  model  over  the  diffusion  approxima¬ 
tion  is  then  discussed. 


H.  THEORY 

The  Boltzmann  equation  for  photon  distribution  function 
7(r,s,0  at  position  r,  direction  s,  and  time  t  from  a  unit 
source  at  position  r0  propagating  along  Sq  at  time  t-  0,  is 
given  by 

d 

~/(r,s,f)  +  cs-  Vr/(r,s,f)  +  c[^(r)  +  ^a(r)]/(r,s,0 

=  c/^(r)J  </s'P(s,s')/(r,s\r)rfs' 

+  S(r-r0)8(s-s0)S(t ),  (2) 

where  c  is  the  speed  of  light  inside  the  medium,  fia  and  fis 
denote  the  position-dependent  absorption  and  scattering  co¬ 
efficients,  and  P(s,s')  is  the  normalized  phase  function  of 
the  light  propagation  in  the  medium.  The  known  phase  func¬ 
tion  is  assumed  to  depend  only  on  the  scattering  angle  ss', 
and  is  then  expandable  in  Legendre  polynomials, 

P(s,s')-(47t)_12  fl/P/(s*s').  (3) 

Equation  (2)  is  nonseparable.  However  the  evolution  in 
direction  space,  F(s,r|s0)exp(-AtflcO  =Jd3r/(r,s,/|r0,So), 
obeys  a  separable  equation  with  the  solution  [14] 

F(s,r|s0)  =  (47r)“12  (ZZ+OexpC-g/OP/is-So).  (4) 

Here  -a,/(2/  + 1)],  especially  g0  =  0  and  g{ 

=  cfifs  where  fi's  is  the  reduced  scattering  coefficient.  The 


formal  solution  to  the  Boltzmann  equation,  Eq.  (1),  is  then 
evaluated  by  (1)  expressing  its  first  <5  function  of  position  r 
as  an  integral  of  exp[/k- (r-c/os(F)<fr')]  over  k  in  the 
Fourier  space,  (2)  making  a  cumulant  expansion  of  the  latter, 
and  (3)  calculating  the  cumulants  in  the  direction  space  with 
use  of  the  exact  Green’s  function  F(s,f|so)  [15]. 

An  arbitrary  order  of  cumulant  solution  can  be  calculated 
[15]  with  higher-order  cumulants  providing  progressively 
more  details  about  the  photon  distribution.  Because  the  pho¬ 
ton  distribution  approaches  a  Gaussian  distribution  when  the 
number  of  the  scattering  events  increases  regardless  of  the 
details  of  the  scattering,  a  second-order  cumulant  solution  is 
sufficient  at  later  times.  At  early  times,  the  photons’  spread  is 
narrow  compared  to  the  resolution  of  the  detector,  hence  the 
detailed  shape  is  less  important  than  the  correct  position  and 
half-width  of  the  beam.  We  emphasize  the  center  of  the  po¬ 
sition  and  half-width  obtained  from  the  second-order  cumu¬ 
lant  solution  is  exact  and  will  not  be  altered  by  higher  order 
cumulant  solutions. 

The  second-order  cumulant  solution  of  the  photon  density 
A^(0)(r,/|r0,s0)  =  /rfs/(0)(r,s,r|r0,s0)  for  an  incident  source 

propagating  along  the  positive  z  axis  (so  =  z)  in  a  uniform 
medium,  is  given  by  [14] 


Ar(0>(r,/|ro,So) 


1  1 


(4iTDzzct) 1/2  4  TrDxxct 


Xexp< 

Xexpj 


(z-;0-/?-)2 
4  D„ct 

(x-x0)2+(y-y0)2 
4  Dxxct 


Xexp  (~/J.act)  (5) 

with  a  moving  center  located  at 

/?.  =  /,[l-exp(-c///,)]  (6) 

and  the  diffusion  coefficients 

Dxx=Dyy 

=f_[l  ,  gzD—  exp(  —  git)]  1—  exp(  —  g2f)’ 
3'151  82i(8i~g2 )  82(81-82)  j’ 

D  =f_[i _ (3gi-g2>n  ~exp(-giQ] 

3/lgi  gi(g  1  82) 


t  2[1  -exp(-g2f)]  3[1  — exp(-git)]2 

82(8 1  82)  2  g] 

For  simplicity,  we  use  the  following  approximation  to  the 
second-order  cumulant  solution  as  the  background  photon 
distribution,  7(0)(r,s,r),  in  an  infinite  uniform  medium  [17], 

/(0)(r,s,t|r0,s0)  =  A^0)(r,t|r0,s0)F(s,t|s0) 

-^D(/)s-V.JV(0>(r,/|ro,So)  (8) 


in  building  the  photon  transport  model  for  image  reconstruc- 
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tion  where  the  time-dependent  diffusion  coefficient  D(t)  is 
taken  to  be  an  average  D(t)  =  (Dxx+Dyy+Dzz)/3  of  the  dif¬ 
fusion  coefficient  ellipsoid.  At  early  times  > 0,  the  first 
term  of  Eq.  (8)  dominates,  and  F(s,f|s0)— ><5(s-s0),  D(t ) 
->c2fV'/9->0,  N^°\r,t\r0,s0)—^S(r—r0—c(t—t0)sQi), 

thus  /(0)(r,s,r|r0,So)  provides  a  correct  picture  of  ballistic 
motion  of  photons  with  speed  c  along  the  incident  direction 

_ I 
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So.  At  later  times,  F(s,r|s0)^(47r)_1,  D(t)—>(3iifs)~l,  Eq. 
(8)  reduces  to  the  photon  distribution  of  the  center-moved 
diffusion  approximation  [8]. 

For  weak  inhomogeneities,  Sjna( r)  and  <5/r'(r),  embed¬ 
ded  in  an  otherwise  uniform  medium,  a  first-order  Bom  ap¬ 
proximation  to  Eq.  (2)  yields  the  change  in  the  photon  dis¬ 
tribution  [17] 


<5/(r,s,Hr0,So)=-^J  dt' j  dr'c8»a(r')N<0\r\t-t%-s)N<0\r\t'\r0,s0)+  ^ J  dt'  j  dr' D(t-t')D(t') 

X[<5/4a(r')+  <5/t'(r')]Vr/A^0)(r',f  — f'|r,  — s)- VrfAf(0)(r',f'|ro,So)+  J  dt' j  dr'D(t-t') 

X[8^a(r')  +  S^'s(r')]exp(-c^'){N^\r',t'\r,-s)s-Vr,N^0\r',t-t'\r0,s0) 

— s0- Vr/^0)(r',r-/'|r,— s)Ar<0>(r',f'|ro,So)}  (9) 


after  neglecting  fast  decaying  terms  involving  exp(-2g,0  for 

1 .  We  should  point  out  that  the  optical  reciprocity  relation 
is  satisfied  by  both  the  photon  density  Eq.  (5)  and  the  photon 
distribution  Eqs.  (8)  and  (9).  At  later  times,  the  term  in  Eq. 
(9)  containing  the  exponential  decay  factor  exp(— cfistf)  can 
be  neglected,  the  change  in  photon  density,  47r<5/(r,s,r),  in 
the  diffusive  limit,  is  reduced  to  that  in  the  diffusion  approxi¬ 
mation  (Eq.  (14)  in  Ref.  [18]). 

The  restriction  of  D(t)  by  taking  an  average  of 
DXX9  Dyy ,  and  Dzz  can  be  relaxed.  The  diffusion  coeffi¬ 
cients  Dxx=Dyy  and  Dzz  can  be  used  instead.  The  only 
change  is  to  replace  all  the  occurrences  of  the  form 
ofI>(OVrriV<0)(r,/|r0,sb)  to 


DJf)  {xdldx +ydldy  )Nw{r,t\ r0  ,s0) 
+  Dzz(/)£o>/<?z'jV(0)(r,t|r0,So) 


in  both  Eqs.  (8)  and  (9). 


A.  Extension  to  planar  geometries 

When  the  scattering  medium  is  bounded,  special  condi¬ 
tions  are  needed  to  set  the  photon  density  at  the  interfaces. 
The  reflection  at  the  interface  reinjects  the  light  into  the  me¬ 
dium.  Using  a  partial  current  technique,  Zhu  et  al  [19] 
showed  that  the  boundary  condition  for  a  semi-infinite  me¬ 
dium  can  be  written  as 

—  0  (10) 

z  =  0 


W<0)- 


dz 


at  the  interface  z  =  0  where 


J2lt  1  ~~Feff 
Z'~3  l+/?eff- 


(ID 


Here  Feff  is  the  effective  reflectivity  at  the  interface  deter¬ 
mined  by  the  Fresnel  reflection  coefficients.  The  extrapola¬ 
tion  length  measures  the  distance  outside  the  medium 
where  the  energy  density  from  the  diffusion  approximation 
vanishes  linearly.  A  recent  study  by  Popescu  et  al  [20]  has 
also  shown  the  dependence  of  the  extrapolation  length  on  the 
scattering  anisotropy. 

The  extrapolated-boundary  condition  has  been  success¬ 
fully  employed  for  planar  geometries  such  as  a  slab  or  a 
semi-infinite  medium  in  diffuse  imaging,  in  which  the  pho¬ 
ton  density  is  set  equal  to  zero  at  an  extrapolated  boundary 
located  a  distance  ze  outside  the  turbid  medium  [8,21,22]. 
The  method  of  images  is  used  to  obtain  the  Green’s  function 
in  such  bounded  media.  The  same  technique  can  be  applied 
here  to  the  Green’s  function  V(0)(r,^|r0,So). 

Keeping  in  mind  that  the  source  approaches  gradually  and 
stops  finally  at  r0+s 0lt  on  average  with  the  increase  of  time, 
the  image  of  the  incident  point  source  at  (-ro^o^o5^) 
propagating  along  the  positive  z  axis  inside  a  semi-infinite 
medium  with  its  interface  at  z  =  0  is  a  negative  one  at 
{xo,yo,-ZQ-2ze-2lt)  propagating  along  the  same  direc¬ 
tion  (Fig.  1).  At  early  times,  both  the  source  and  its  image 
have  not  arrived  at  the  extrapolated  boundary  and  their  con¬ 
tributions  at  the  extrapolated  boundary  can  be  neglected. 
When  the  time  increases,  the  contributions  at  the  extrapo¬ 
lated  boundary  from  both  the  source  and  image  tend  to  can¬ 
cel  each  other  as  both  approach  their  final  stops  (shadow 
spots  in  Fig.  1).  The  shadow  spots  just  represent  the  positions 
of  the  source  and  its  image  in  the  center-moved  diffusion 
approximation.  The  Green’s  function  of  a  semi-infinite  me¬ 
dium  given  by 
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FIG.  1.  The  incident  source  at  position  (jc0,y0,z0>0)  and  its 
image  source  at  (x0,y0,-zo~2ze—2ll)  propagating  along  the 
positive  z  axis  in  a  semi-infinite  medium  (z>0)  with  its  interface  at 
z  =  0.  The  source  and  its  image  move  from  their  original  positions 
(dark  spots)  to  their  final  stops  (shadow  spots)  at  later  times. 

^semi(r’^lrO’So)  >yo  »Sq) 

-N<0\r,t\x0,y0,-z0-2ze-2ll  ,Sq) 

(12) 

thus  approximately  satisfies  the  extrapolated-boundary  con¬ 
dition. 

The  same  procedure  can  be  easily  applied  to  a  slab  with 
its  extrapolated  boundaries  at  z  =  0  and  z  =  L.  The  images  of 
an  incident  source  at  U0^o^o)  with  propagating 

along  positive  or  negative  z  axis  (sz=  ±  1 )  are  a  set  of  posi¬ 
tives  images  at  (xQiy0,Zo  +  2nL)  and  a  set  of  negative  ones 
at  (x0,y0,-z-2nL-2szlt),  all  propagating  along  the  same 
direction  as  the  source  (  — oo<w<oo  is  integer). 

B.  Comparison  with  the  Monte  Carlo  simulation 

We  will  compare  the  photon  densities  computed  by  the 
diffusion  approximation  (DA),  the  cumulant  approximation 
(CA)  Eqs.  (5)  and  (8),  and  the  Monte  Carlo  method  (MC)  for 
an  incident  collimated  pulse  first  in  an  infinite  space  and  then 
in  a  semi-infinite  space.  In  DA,  the  incident  photons  are  as¬ 
sumed  initially  scattered  isotropically  at  a  depth  of  one  trans¬ 
port  mean  free  path  into  the  medium  as  used  by  Patterson 
et  al.  [8].  No  such  adjustments  are  performed  in  CA.  The 
Monte  Carlo  code  is  adapted  from  Prahl  et  ai  [23]  and  Wang 
et  al  [24].  Photons  are  launched  one  by  one  into  the  me¬ 
dium.  Each  photon  (regarded  as  a  packet)  starts  from  the 
origin  of  the  coordinate  system  and  the  first  scattering  event 
takes  place  along  the  positive  z  axis.  The  step  size  (distance 
between  consecutive  scattering  events)  is  sampled  from  an 
exponential  distribution  characterized  by  the  total  attenuation 
+  following  Beer’s  law.  After  each  propagation 
step,  the  photon  packet  is  split  into  two  parts — a  fraction 
(Ma  / Mr)  is  absorbed  and  the  rest  scattered.  The  new  propa¬ 
gation  direction  after  scattering  (three  directional  cosines)  is 
sampled  by  assuming  a  Henyey-Greenstein  phase  function 
[25].  The  effect  of  internal  reflection  is  included  in  this  code. 
The  technique  of  roulette  [26]  is  used  to  terminate  a  photon 


FIG.  2.  The  center  position  and  the  half-width  of  the  photon 
cloud  inside  a  uniform  infinite  absorptionless  medium  with  aniso¬ 
tropy  equal  to  0.9. 


packet  to  improve  the  efficiency  of  the  calculation  without 
introducing  a  bias.  The  results  in  the  following  paragraphs 
have  been  scaled  to  use  the  transport  mean  free  path  /,  as  the 
unit  of  the  length  and  the  flight  time  for  one  transport  mean 
free  path  in  the  medium  lt/c  as  the  unit  of  the  time.  The 
source  is  incident  along  the  positive  z  axis  at  the  origin  in 
space  and  time.  5  X 106  photons  are  used  in  one  run  of  the 
Monte  Carlo  simulation. 

The  first-  and  second-order  cumulants  (the  center  position 
and  the  half- width  of  the  “photon  cloud”)  of  our  cumulant 
solution  Eq.  (5)  is 

(^(f))==(y(0)=0,  (z(f))=/,[l  -exp(-cr//,)], 


]2Dxx(t)ct, 


where  ()  means  an  ensemble  average  of  photon  positions  at  a 
specified  time.  This  theoretical  prediction  can  be  easily  veri¬ 
fied  by  a  Monte  Carlo  simulation.  Figure  (2)  shows  the  first 
two  cumulants  of  photons  for  an  incident  pulse  along  the  z 
axis  at  time  zero  into  an  infinite  medium  with  anisotropy  0.9. 
A  perfect  agreement  on  the  center  position  (the  first-order 
cumulant)  and  the  half-width  (the  second  order  cumulant)  of 
the  photon  distribution  between  our  theoretical  result  and  the 
Monte  Carlo  simulation  is  obtained.  The  half-widths  along 
xyz  directions  are  very  close;  the  value  along  z  direction  is  a 
bit  larger  than  that  along  the  xy  direction  as  predicted  by  Eq 
(13). 

Figures  3(a)-3(c)  shows  the  photon  density  at  positions 
(0,0,3/,),  (0,0,6/,),  and  (0,0,10/,)  computed  by  all  three  dif¬ 
ferent  methods  for  the  same  infinite  medium.  At  a  distance  of 
3/,,  the  time  profile  of  photon  density  from  the  cumulant 
approximation  agrees  much  better  to  the  Monte  Carlo  result 
than  DA  by  providing  a  correct  peak  position  of  photon  den¬ 
sity.  Some  amount  of  photons  arriving  faster  than  the  speed 
of  light  still  exist  in  this  second-order  cumulant  calculation. 
However,  this  is  already  a  big  improvement  compared  to 
DA.  The  result  from  CA  can  be  further  improved  when 
higher-order  cumulants  are  used.  At  a  larger  distance,  all  the 


066609^ 


PHOTON  MIGRATION  IN  TURBID  MEDIA  USING  A  . . . 


FIG.  3.  Photon  density  at  positions  (a)  (0,0,3/,),  (b)  (0,0,6/,), 
and  (c)  (0,0,10/,)  vs  time  normalized  to  a  unit  source  in  an  infinite 
medium.  The  source  is  incident  along  the  positive  z  axis  at  the 
origin  of  the  coordinate  system  and  at  time  zero.  The  three  curves 
are  computed  by  the  diffusion  approximation  (DA),  the  cumulant 
approximation  (CA),  and  the  Monte  Carlo  method  (MC),  respec¬ 
tively. 
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FIG.  4.  Photon  density  at  positions  (a)  (0,0,3/,),  (b)  (0,0,6/,), 
and  (c)  (0,0,10/,)  vs  time  normalized  to  a  unit  source  in  a  semi¬ 
infinite  medium.  The  source  is  incident  normal  to  the  surface  of  the 
medium  and  along  the  positive  z  axis  at  the  origin  of  the  coordinate 
system  and  at  time  zero. 
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FIG.  5.  The  backscattered  photon  intensity  I(r,-z,t)  at  positions  (a)  (0,/„0)  and  (b)  (0,2/„0)  on  the  boundary  of  a  semi-infinite  medium 
vs  tune  normalized  to  a  unit  source  in  a  semi-infinite  medium.  The  source  is  incident  normal  to  the  surface  of  the  medium  and  along  the 
positive  z  axis  at  the  origin  of  the  coordinate  system  and  at  time  zero. 


three  methods  begin  to  agree  with  each  other  pretty  well  and 
the  cumulant  approximation  is  better  than  the  diffusion  ap¬ 
proximation.  The  difference  between  results  from  DA,  CA, 
and  MC  is  negligible  when  the  distance  is  10/,  or  larger. 

The  calculations  using  DA,  CA,  and  MC  are  performed 
again  for  a  semi-infinite  medium  with  its  boundary  at  z  =  0, 
whose  optical  parameters  are  taken  to  be  the  same  as  the 
above  infinite  medium.  The  source  is  incident  at  the  origin  of 
the  coordinate  system  and  along  the  positive  z  axis  (normal 
to  the  surface)  at  time  zero.  The  effective  reflectivity  is  taken 
to  be  zero  and  an  extrapolation  length  z*=0.7/,  is  used  in 
both  DA  and  CA  calculations.  Figures  4(a)-4(c)  shows  the 
corresponding  results  for  this  example.  Again  CA  shows  a 
much  better  agreement  to  the  MC  than  DA.  Compared  with 
Fig.  3  for  the  infinite  case,  the  tail  of  the  profiles  in  Figs. 
4(a)-4(c)  for  the  semi-infinite  medium  is  lower  due  to  the 
presence  of  an  extra  negative  image  source  coming  from  the 
boundary  condition. 

As  a  final  example,  the  backscattered  photon  intensity 
7(0)(r,~ £,/)  at  positions  (0,/„0)  and  (0,2/„0)  on  the  bound¬ 
ary  of  the  above  semi-infinite  medium  is  calculated  with  use 
of  the  three  different  methods  (see  Fig.  5).  In  DA,  photons 
diffuse  from  the  adjusted  source  position  (0,0,/,)  with  the 
constant  diffusion  coefficient  D  =  /,/3.  In  CA,  back-scattered 
photons  arrive  later  because  the  center  of  photons  moves 


forward  along  the  positive  z  direction  and  diffuse  from  the 
moving  center  with  a  gradually  increasing  diffusion  coeffi¬ 
cient  from  0  to  /,/ 3.  CA  agrees  well  with  the  Monte  Carlo 
simulation. 

C.  Weight  function  for  image  reconstruction 

The  response  (the  change  of  the  scattered  wave  field)  to  a 
unit  absorption  or  scattering  inhomogeneity  is  usually  called 
the  weight  function  or  the  Jacobian  in  medical  tomography 
literature.  This  quantity  plays  a  central  role  in  image  recon¬ 
struction  regardless  of  which  method  is  used  to  obtain  the 
inhomogeneity  distribution  in  a  medium.  Let  us  rewrite  Eq. 
(9)  in  the  following  form: 

4TrS/(r,s,tlz0,s0)=-cf  dr'  S/ia(r')wa(r,s,r0,So,t;r') 


Xw’f(r,s,r0,s0,f;r')  (14) 

with  the  absorption  and  scattering  weight  functions  defined 
as 


wa(r,s,r0,So,f;r')=  Jo^'^0>(r',r-r'|r,-s)^<0)(r'(/'|r0,s0)-W,(r,s,r0,So,r;r')/(3^2) 

ws(r,s,r0,So,r;r')  ft  rt 

- — - =  U/'D(f-r')£>(/')Vr^(0>(r',f-r'|r,-s)-Vr,/V<0)(r',/'|r0,s0)+  dt'D(t') 

VPs  Jo 

Xexp[— £•/*,'(?  — /')]JV<0)(rV  —  /'|r,  —  s)s-Vr/N(0)(r',f'|ro,so)—  ['dt'D(t-t') 

Jo 

Xexp(  — c/z'f')so- Vr-W(0>(r',f— /'|r,  — s)Af(0)(r',f'|ro,So), 


05) 
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respectively. 

As  A^0)(r,f|r0,So)— ><5(r— r0— cfso)  when  t->  0,  special  attentions  must  be  paid  when  a  numerical  integration  is  carried  out 
for  Eq.  (15).  The  range  of  integration  is  divided  into  three  areas:  (0,A),  (A, I- A),  and  (t-A,t)  where  />A>0.  The  end 
corrections  from  the  integration  over  (0,A)  and  (f- A,f)  to  the  weight  functions  integrated  over  (A,f- A)  range  are  approxi¬ 
mately  given  by 


efl(r,s,ro,so,t;r')=-^(r'(t|r-s)5(x'-xo)%'-yo)mA-aH(|)+^(r',t|ro,So)5(x'-^)<5(y'-y)H(A-17)H(77) 


-^(r,s,r0,So,t;r')/(3/A'2), 


^(r.s.ro.so.fjr')  _  D(t)  (  ^  dDjt) 


9/t. 


r  2 


x  l 


mt) 


cdt 


cdt 


t=v 


t=£‘ 


s0-VrMr\t\r-s)S(x'-x0)S(y,-y0)H(A-Z)H(Z)  + 


D(t) 


s- VrrN(rf 9t\r0,So)S(x' -x)S(y’ -y)H(A-  tj)H{ tj ), 


(16) 


where  r'  =  r0+c£s0  =  r-  77CS  is  the  position  of  the  inhomo¬ 
geneity,  s  and  So  are  assumed  to  be  along  the  positive  or 
negative  z  axis,  and  H  is  the  Heaviside  function.  In  the  fol¬ 
lowing  calculations,  the  refractive  index  of  the  medium  is 
assumed  to  be  1.33;  the  absorption  and  scattering  coefficients 
of  the  medium  are  assumed  to  be  0.003  mm  1  and 
10  mm”1,  respectively,  and  the  scattering  anisotropy  0.9, 
providing  a  transport  mean  free  path  /,=  1  mm.  The  offset  A 
is  taken  to  be  0.1  ps  when  more  than  99%  of  the  photon 
packet  still  concentrates  within  a  cubic  volume  (0.01/, )3. 
The  weight  function  from  a  cubic  of  volume  (0.01/, )3  is 
calculated  using  the  DCUHRE  algorithm  [27]. 

The  weight  functions  for  a  semi-infinite  medium  are 
shown  in  Fig.  6  for  absorption  and  scattering  inhomogene¬ 
ities.  The  backscattered  photons  (propagating  along  negative 
z  axis)  are  detected  by  a  detector  placed  at  a  position 
(0,2/„0),  off  two  transport  mean  free  path  from  the  source. 
Figures  6(a)  and  6(c)  show  the  response  to  an  inhomogeneity 
at  (0,0, z)  positions  which  is  in  the  propagation  direction  of 
the  source  at  delay  50  ps  and  500  ps.  The  CA  shows  a  much 
stronger  response  from  the  inhomogeneity  in  the  propagation 
direction  of  and  close  to  the  source  than  the  diffusion  ap¬ 
proximation.  Both  absorption  and  scattering  weight  functions 
from  CA  reveal  a  peak  at  about  0.03/, .  This  peak  originates 
from  the  initial  ballistic  motion  of  the  incident  photon.  In  a 
short  time  after  the  photon  is  launched  (f— > 0),  the  photon 
packet  will  be  positioned  at  z*  =  ct  with  a  spread  of  half¬ 
width  ~  yj4D(t)ct=2ctylct/(3Jf),  hence  the  presence  of  an 
absorption  or  scattering  inhomogeneity  at  position  (0,0, z 
<z*),  sitting  in  the  ballistic  path  of  the  photon,  will  signifi¬ 
cantly  reduce  the  number  of  backscattered  photons  received 
by  the  detector  [wa>0  and  ws<0  in  Eq.  (14)].  In  Figs.  6(b) 
and  6(d)  where  the  inhomogeneity  is  placed  at  (0 ,/,,z)  mm 
positions,  not  sitting  in  the  ballistic  path  of  the  photon,  this 
peak  is  gone. 

The  diffusion  approximation  is  invalid  when  the  inhomo¬ 
geneity  is  too  close  to  the  source  or  the  detector.  Neverthe¬ 
less,  the  weight  function  from  DA  is  plotted  over  the  full  z 


I 

range  for  comparison.  A  peak  is  found  in  the  absorption 
weight  function  [Fig.  6(a)]  and  a  crossing  of  zero  in  the 
scattering  weight  function  [Fig.  6(c)]  at  z  =  /, ,  because  of  the 
artificial  adjustment  of  the  source  position  to  one  transport 
mean  free  path  into  the  medium  and  the  singularity  of  the 
Green’s  function  in  DA  when  the  inhomogeneity  and  the 
source  overlap. 

A  larger  disagreement  between  C A  and  DA  is  observed  in 
the  scattering  weight  function  than  in  the  absorption  weight 
function.  The  absorption  weight  functions  from  CA  and  DA 
agree  with  each  other  relatively  well  except  for  a  region  of 
depth  of  /,  near  the  surface  when  the  inhomogeneity  is  in  the 
propagation  direction  of  the  source,  or  in  the  field  of  view  of 
the  detector.  The  scattering  weight  functions  from  CA  and 
DA  disagree  significantly  within  the  region  of  depth  of  at 
least  2/,  close  to  the  surface.  The  deepest  position  that  can  be 
detected  by  the  detector  at  time  t  is  roughly  ct! 2.  This  con¬ 
dition  is  better  observed  by  CA  because  CA  shows  a  faster 
decay  rate  of  both  the  absorption  and  scattering  weight  func¬ 
tions  with  the  increase  of  the  depth  [Figs.  6(a)-6(d)]. 

The  absorption  and  scattering  weight  functions  for  a  slab 
is  shown  in  Fig.  7.  The  slab  has  the  same  optical  parameter 
as  the  semi-infinite  medium.  The  thickness  of  the  slab  is  d 
=  30/,.  The  source  is  at  the  origin  (0,0,0).  The  detector  is 
placed  on  the  opposite  side  of  the  slab,  (0,0,30/,),  in  the 
propagation  direction  of  the  source.  The  weight  functions  by 
the  cumulant  approximation  and  the  diffusion  approximation 
are  strictly  symmetric  about  the  plane  z  =  d/2.  The  agreement 
between  CA  and  DA  for  the  absorption  weight  function  is 
better  than  for  the  scattering  weight  function.  Both  CA  and 
DA  produce  close  results  for  the  inhomogeneity  not  located 
near  the  boundary.  When  the  inhomogeneity  is  placed  along 
the  line  (0,0,z),  in  the  propagation  direction  of  the  source 
and  in  the  field  of  view  of  the  detector,  two  peaks  at  about 
0.03/,  and  d-0.03/,  appear  in  CA;  two  peaks  in  the  absorp¬ 
tion  weight  function  [Fig.  7(a)]  and  two  crossings  of  zero  in 
the  scattering  weight  function  [Fig.  7(c)]  appear  in  DA  at  /, 
and  d— /,. 
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(a) 


(b) 


z/l, 


<C)  (d) 

FIG.  6.  Weight  functions  for  a  semi-infinite  medium  where  the  inhomogeneity  is  (a)  absorption  and  (c)  scattering  at  (00z)  in  the 
propagation  direction  of  the  source;  (b)  absorption  and  (d)  scattering  at  (0,/,  ,*),  off  by  one  transport  mean  free  path.  Profiles  at  two  delay 

times  f-50  ps  and  t - 500  ps  are  plotted  for  both  the  cumulant  approximation  (CA)  and  the  diffusion  approximation  (DA).  The  insets  replot 
the  weight  functions  in  a  logarithm  scale.  F 


m.  DISCUSSION 

We  attribute  the  formation  of  a  peak  very  close  to  the 
surface  but  not  on  the  surface  (about  0.031 1  into  the  medium) 
of  the  absorption  and  scattering  weight  function  in  the  cumu¬ 
lant  approximation  to  the  initial  ballistic  motion  of  the  inci¬ 
dent  photon.  The  photon  penetrates  into  the  medium  with  an 
initial  speed  of  c  and  with  its  center  approaching  and  stop¬ 
ping  at  one  transport  mean  free  path  into  the  medium.  Hence 
the  effect  is  only  significant  when  the  inhomogeneity  is  in 
the  propagation  direction  of  the  source  or  in  the  field  of  view 
of  the  detector,  and  the  peak  response  shifts  away  from  the 
surface  of  the  medium. 

The  diffusion  approximation  requires  one  to  adjust  the 
position  of  the  source  to  compensate  for  the  initial  ballistic 
motion  of  the  photon  [8,28,29].  From  our  more  accurate  re¬ 
sult  Eq.  (9),  the  source  and  the  detector  terms  appear  in  a 
form  of  W<0)(r  ,r|r,  — s)  and  A^0)(r',7jr0,So),  respectively. 
The  source  and  the  detector  approach  gradually  and  stop  at 
r0+/^o  311(1  r-lfi,  respectively,  with  the  increase  of  time 


where  Sq  and  s  are  the  propagating  directions  of  the  incident 
and  outgoing  photon.  The  positioning  of  both  the  source  and 
the  detector  for  one  transport  mean  free  path  into  the  me¬ 
dium  is  hence  mandatory  if  the  diffusion  approximation  is 
used.  The  curves  for  DA  in  Figs.  6  and  7  are  calculated  using 
this  adjustment.  The  DA  will  deviate  from  the  CA  signifi¬ 
cantly  over  the  full  range  of  depth  if  the  adjustment  on  the 
position  of  the  source  or  the  detector  is  not  performed. 

The  diffusion  approximation  for  image  reconstruction 
substantially  underestimates  the  contribution  to  the  emission 
measurement  from  the  inhomogeneity  in  the  propagation  di¬ 
rection  of  and  close  to  the  source,  or  in  the  field  of  view  and 
close  to  the  detector.  This  error  may  distort  the  signal  from 
the  inhomogeneity  inside  the  medium  because  the  weight 
function  near  surface  is  usually  much  larger  than  that  inside, 
and  may  lead  to  a  failure  in  image  reconstruction.  The  high 
response  from  the  region  near  surface  is  not  desirable  when 
the  inhomogeneity  inside  the  medium  is  to  be  detected  in  the 
transmission  or  backscattering  measurements.  The  cancella- 
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FIG.  7.  Weight  function  for  a  slab  where  the  inhomogeneity  is  (a)  absorption  and  (c)  scattering  at  10,0,-).  in  the  propagation  direction  of 
the  source,  (b)  absorption  and  (d)  scattering  at  (0 ,/,,z),  off  by  one  transport  mean  free  path.  Weight  functions  calculated  from  the  cumulant 
(CA)  and  diffusion  (DA)  approximations  are  plotted  for  time  delays  of  r=300  ps  and  t=  1500  ps.  The  insets  replot  the  weight  functions  in 
a  logarithm  scale. 


tion  between  multiple  measurements  using  nearby  wave¬ 
lengths  may  help  reduce  this  effect. 

Other  attempts  to  obtain  a  better  approximation  to  radia¬ 
tive  transfer  in  turbid  media  were  made  by  various  authors 
such  as  Ishimaru’s  diffusion  approximation  [30,31],  the  te¬ 
legrapher  equation  of  Durian  et  al.  based  on  the  two  stream 
theory  [32],  Gershenson’s  time-dependent  equation  in  the 
diffusion  limit  using  a  higher-order  angular  expansion  [6], 
and  non-Euclidean  diffusion  equation  of  Polishchuk  et  al. 
[33].  The  advantage  of  this  cumulant  approximation  is  that  it 
provides  a  clear  picture  of  photon  migration  for  an  incident 
collimated  beam  from  early  to  later  times  and  that  it  gives 
the  exact  center  and  the  exact  spread  of  the  photon  cloud  at 
both  early  and  later  times  by  only  using  a  second-order  cu¬ 
mulant  approximation. 

In  conclusion,  we  have  presented  a  cumulant  approxim¬ 
ation  to  radiative  transfer,  which  provides  an  analytical  tool 
to  describe  photon  migration  at  both  early  and  later  times — 


from  the  initial  ballistic  motion  till  the  final  diffuse  regime. 
To  a  second-order  cumulant,  the  solution  agrees  with  the 
Monte  Carlo  simulation  at  later  times  and  provides  a  correct 
peak  position  in  time  for  photon  arrivals  at  early  times,  in 
both  an  infinite  medium  and  a  bounded  medium  with  a  pla¬ 
nar  geometry.  The  initial  ballistic  motion  of  photon  produces 
a  strong  peak  in  the  response  from  absorption  and/or  scatter¬ 
ing  inhomogeneities,  which  are  in  the  propagation  direction 
of  and  close  to  the  source,  or  in  the  field  of  view  of  and  close 
to  the  detector,  at  both  early  and  later  times. 
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ABSTRACT 

Spectroscopic  and  time-sliced  two-dimensional  (2-D)  transillumination  imaging  methods  were  used  to  investigate  ex 

vivo  tumor  and  normal  tissues  of  human  breast  and  parotid  gland.  The  experimental  arrangement  for  time-sliced  optical 

imaging  uses  120-fs,  1  kHz  repetition-rate,  800-nm  light  pulses  from  a  Tirsapphire  laser  system  for  sample  illumination 

and  a  charge  coupled  device  (CCD)  camera  coupled  to  a  gated  image  intensifier  for  recording  2-D  images.  The 

spectroscopic  imaging  arrangement  uses  1210-1325  nm  tunable  output  of  a  Cr:  forsterite  laser  for  sample  illumination,  a 

Fourier  space  gate  to  discriminate  against  multiple-scattered  light,  and  a  near-infrared  (NIR)  area  camera  to  record  2-D 

images.  Images  recorded  with  earlier  temporal  slices  of  transmitted  light  highlighted  tumors,  while  those  recorded  with 

later  slices  accentuated  normal  tissues.  When  light  was  tuned  closer  to  the  1203  nm  absorption  resonance  of  adipose 

tissues,  a  marked  enhancement  in  contrast  between  the  images  of  adipose  and  fibrous  tissues  was  observed.  A  similar 

wavelength-dependent  difference  between  normal  and  cancerous  tissues  was  observed.  These  results  correlate  well  with 

pathology  and  nuclear  magnetic  resonance  based  analyses  of  the  samples.  This  work  demonstrates  the  advantages  of 

time-resolved  spectroscopic  approach  for  imaging  tumors  in  body  organs.  ' 

Keywords :  Time-resolved  optical  imaging,  near-infrared  spectroscopic  imaging,  breast  cancer,  nuclear  magnetic 

resonance,  noninvasive  cancer  detection,  parotid  gland,  Warthin’s  tumor 

1.  INTRODUCTION 

The  need  for  noninvasive  imaging  modalities  for  early  detection  and  diagnosis  of  cancer  drives  the  recent  research  * 

interest  in  developing  NIR  light  based  imaging  techniques.1'4  Many  researchers  use  measurement  of  amplitude  and 
phase  of  photon  density  wave  launched  in  a  tissue  specimen  upon  irradiation  with  frequency-modulated  NIR  light  to 
obtain  images  of  the  target  specimen.4  We  have  developed  a  fast,  time-resolved  imaging  approach5"7  because  it  is 
intuitive  simple  and  can  provide  simultaneous  information  at  a  multitude  of  frequencies  that  is  useful  for  three- 
dimensional  image  reconstruction.  We  complement  the  time-sliced  imaging  at  800  nm  with  continuous  wave  (CW) 
spectroscopic  imaging  that  uses  NIR  light  in  the  therapeutic  wavelength  range  of  1200-1350  nm  to  realize  the  diagnostic 
potential  of  optical  imaging  approach.8  j 

An  important  step  in  the  development  of  a  noninvasive  optical  imaging  approach  is  testing  its  efficiency  on  ex-vivo 
tissue  specimens  of  the  target  organ,  which  in  turn  provides  key  information  about  light  transport  and  optical 
characteristics  of  the  types  of  tissues  under  investigation.  In  this  article  we  present  result  of  measurements  on  human 
breast  and  parotid  gland  tissues. 

2.  EXPERIMENTAL  ARRANGEMENT 

Our  experimental  arrangement  for  time-sliced  imaging,  detailed  elsewhere,5'7  made  use  of  approximately  120-fs 
duration,  1-kHz  repetition-rate,  800-nm  light  pulses  form  a  Ti-sapphire  laser  and  amplifier  system9  for  sample 
illumination,  and  an  ultrafast  gated  intensified  camera  system  (UGICS)  for  recording  2-D  images.7  The  UGICS  system 
provides  a  80-ps  time  gate  whose  position  can  be  varied  over  a  20  ns  range.7  The  arrangement  for  NIR  spectroscopic 
imaging  used  die  output  of  a  Cr4+:  forsterite  laser  to  illuminate  the  sample.  A  Fourier  space  gate10  selected  out  a  fraction 
of  the  less  scattered  image-bearing  photons,  and  a  50  mm  focal-length  camera  lens  directed  them  to  the  320  x  240  pixels 
imaging  element  of  an  InGaAs  NIR  area  camera.  A  1 0  mm  x  20  mm  x  5  mm  thick  human  breast  tissue  specimen 
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composed  of  normal  adipose  (NA),  and  cancerous  (C)  tissues  was  used  in  the  experiment.  The  sample  obtained  from  the 
left  breast  of  a  63-year  old  female  following  a  modified  radical  mastectomy  and  lymph  node  dissection  at  the  Memorial 
Sloan-Kettering  Cancer  Center  (MSKCC)  was  provided  to  us  under  an  Institutional  Review  Board  (IRB)  approval  from 
CCNY.  The  clinical  diagnosis  revealed  an  invasive  ductal  carcinoma  (IDC)  with  poorly  differentiated  histologic  grade 
III/III  and  nuclear  grade  III/III.  The  specimen  was  placed  in  a  high  grade  Pyrex  glass  (suitable  for  both  optical  and 
NMR  studies)  container  that  slightly  compressed  and  retained  the  tissue  axial  orientation.  Following  the  optical 
imaging,  the  undisturbed  samnle  cell  was  transferred  to  MSKCC  and  weighted  Tx  and  T2  relaxation  time  NMR  images  of 
the  sample  were  recorded.11*1* 

The  20  mm  x  20  mm  x  10  mm  excised  parotid  gland  tissue  specimen  was  obtained  from  the  New  York  Eye  and  Ear 
Infirmary  (NYEEI),  under  an  IRB  approval  from  CCNY.  The  sample  was  placed  between  two  glass  slides  and  mildly 
compressed  to  maintain  a  uniform  thickness.  The  specimen  came  from  the  left  parotid  gland  of  a  37-year-old  male 
patient.  The  clinical  diagnosis  revealed  a  Warthin’s  tumor13  (WT)  of  the  left  parotid  gland  (PG). 

3.  EXPERIMENTAL  RESULTS 


3.1  Time-sliced  imaging 


Time-sliced  2-D  trans illumination  images  of  the  breast  tissue 
sample  for  gate  positions  of  25-ps  and  225-ps  are  displayed  in  the 
upper  frames  of  Figs.  1(a)  and  1(b),  respectively.  The  zero 
position  was  taken  to  be  the  time  of  arrival  of  a  light  pulse  through 
a  5-mm  thick  quartz  cell  filled  with  water.  The  corresponding 
lower  frames  present  the  spatial  intensity  profiles  of  the  respective 
images  integrated  over  the  same  horizontal  area  in  both  the  images. 
The  salient  feature  of  the  images  is  the  difference  in  time- 
dependent  brightness  of  the  normal  tissue  (N)  and  cancer  (C) 
regions  in  the  sample.  In  the  25-ps  image  of  Fig.  1(a)  the  cancer 
region  (C)  is  prominent.  The  corresponding  horizontal  spatial 
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Fig.  1  Time-sliced  images  of  the  breast  tissue  specimen  %  120  a  ;  60  ; 

at  the  gate  position  of  (a)  25  ps,  and  (b)  225  ps.  J  ^  j M  l  jr\ 

intensity  profiles  exhibit  peak  in  the  cancer  (C)  *  J  \ ;  30.  A  \ 

region.  At  this  early  time,  markedly  more  light  is  |  40  j  WT  I;  NPG  WT  \ 

transmitted  through  the  tumor  than  through  the  2  n  / _ a  j  j NPG  \ 

normal  tissue.  With  time  the  normal  tissue  image  U  100  ^  200  °  100  150  200 

gain  in  brightness  as  shown  in  the  225-ps  image  of  Pixels  Pixels 

Fig.  1(b).  The  corresponding  spatial  intensity  Fig.  2  Time-sliced  images  of  the  parotid  gland  tissue  specimen 

profile  peaks  in  the  normal  (N)  position  indicating  at  the  gate  position  of  (a)  50  ps,  and  (b)  250 ps. 

higher  light  transmission  through  the  normal  region 

than  the  tumor.  This  result  indicates  that  light  transits  faster  through  tumor  than  normal  tissue.  Lower  scattering  or/and 
higher  absorption  of  light  by  the  tumor  may  account  for  the  observed  temporal  behavior.  Since,  there  is  no  known 
absorption  of  800-nm  light  by  the  breast  tissue,  we  attribute  these  time-dependent  differences  in  transmitted  light 
intensity  to  higher  scattering  of  light  by  normal  tissue  than  the  cancerous  tissue. 


Similar  time-dependent  differences  in  transmitted  light  intensity  and  image  characteristics  are  observed  between  the 
Warthin’s  tumor  and  the  normal  regions  of  the  parotid  gland  tissue  sample,  as  shown  in  Figs.  2(a)  and  2(b). 

3.2  Spectroscopic  imaging 

Spectroscopic  images  of  the  breast  tissue  sample  recorded  using  light  of  wavelength  1225  nm,  and  1300  nm  appear  in 
the  upper  frames  of  Figs.  3(a)  and  3(b),  respectively.  Corresponding  spatial  profiles  integrated  over  the  same  horizontal 
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Fig.  4  Spectroscopic  images  of  the  parotid  gland  specimen 
at  the  wavelength  of  (a)  1 225  nm,  and  (b)  1285  nm. 

area  in  all  images  appear  in  the  respective  lower  frames.  The 

s«  mi  w  50  100  150  salient  feature  of  the  images  is  the  higher  brightness  of  the 

Pixels  Pixels  —  - 

Fig.  3  Spectroscopic  images  of  the  breast  tissue  sample 
at  the  wavelength  of  (a)  1225  nm,  and  (b)  1300  nm. 

cancer  (C)  region  compared  to  the  normal  (N)  region  for 
both  the  wavelengths.  This  indicates  higher  overall  light 
transmission  through  the  tumor  than  through  normal 
tissue.  A  more  interesting  parameter  is  the  average 
value  of  the  tumor-to-normal  intensity  ratio  RtnM  as  a 
function  of  wavelength.  The  value  of  /?tn  is 
approximately  2.2  at  1225  nm  and  1.5  at  1300  nm,  a  Fi8-  5  NMR  ima8es  ofTi  relaxation  time  of  the  5-mm  thick  breast 

significant  difference.  tissue  samPle’  (a)  ^rst  2.5-mm  axial  section,  and  (b)  second 

2.5-mm  axial  section. 

Similar  behavior  is  observed  in  the  spectroscopic  images  of  the  Warthin’s  tumor  (WT)  and  normal  parotid  gland  (NPG) 
tissues  displayed  in  Figs.  4(a)  and  4(b).  The  value  of Rjs  is  approximately  2.9  at  1225  nm  and  1 .7  at  1285  nm,  again  a 
significant  difference. 

33  Nuclear  magnetic  resonance  imaging 

NMR  Ti- weighted  images  of  the  breast  tissue  specimen  are  shown  in  Fig.  5(a)  and 
5(b).  Each  figure  represents  an  axial  planar  image  of  a  2.5-mm  thick  region  of  the 
sample.  Repetition  time  (TR)  for  obtaining  these  images  was  500-ms  which  is 
greater  than  300-ms,  the  lower  limit  required  for  a  good  contrast/noise  ratio 
(CNR).[13]  Echo  time  (TE)  was  16-ms  much  shorter  than  80-ms  that  also  implies  a 
fairly  good  effective-contrast. [13]  The  T^-weighted  images  showed  no  further 
improvements  in  image  contrast. 

Fig.  6  Histological  map  of  the  breast  For  even  ftlrther  validation  of  the  optical  imaging  results,  we  obtained  histological 
tissue  sample  Siowingnormal  adipose  micrographs  of  both  the  breast  and  the  parotid  gland  samples.  Figure  6  shows  a 

(NA),  and  cancerous  (C)  region.  histological  micrograph  of  the  breast  tissue  sample.  The  upper  left  part  that  appear 

much  denser  than  the  rest  of  the  micrograph  has  an  abundance  of  cancer  cells  (C), 
the  lower  left  part  has  extra-cellular  collagen  matrix  with  some  interspersed  cells,  and  the  much  lighter  upper  right  part  is 
rich  in  normal  adipose  (NA)  tissues.  The  histology  results  are  consistent  with  the  optical  imaging  and  NMR  results. 
Similar  validation  of  optical  imaging  results  for  the  parotid  gland  tissue  specimen  was  also  obtained  from  histological 
analysis. 

4.  DISCUSSION 


The  results  of  both  the  time-sliced  imaging  and  spectroscopic  imaging  experiments  show  the  ability  of  the  optical 
imaging  methods  to  select  between  the  tumor  and  normal  tissues.  Images  recorded  with  the  earlier  arriving  light 
highlighted  the  tumor,  while  those  recorded  with  later  arriving  light  accentuated  the  normal  region.  The  normal  upper 
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right  part  of  the  breast  specimen  was  mostly  adipose  tissue.  However,  we  have  observed  similar  time-dependent 
difference  in  light  transmission  between  cancer  and  normal  fibrous  breast  tissues  as  well.[5,14]  Our  observation  is  that 
for  the  same  thickness  of  tissues,  light  arrives  earliest  through  the  cancerous  tissue,  next  through  normal  fibrous  tissues, 
and  latest  through  the  adipose  tissues  in  a  breast  tissue  specimen.  Similar,  time-dependent  transmission  was  observed 
for  Warthin’s  tumor  and  normal  parotid  gland  tissues. 

Spectroscopic  imaging  results  are  encouraging  and  may  be  useful  in  obtaining  diagnostic  information.  In  particular,  the 
wavelength-dependent  tumor-to-normal  transmitted  intensity  ratio  Rm(X)  has  the  potential  to  be  a  useful  parameter  in 
tumor  identification.  Extensive  data  with  different  types  of  cancer  will  be  needed  to  establish  this  potential. 

Validation  of  optical  imaging  results  is  an  important  consideration.  While  pathology  is  ideal  for  validation  of  ex  vivo 
measurements,  an  alternative  is  needed  for  in  vivo  applications.  Magnetic  resonance  imaging  carried  out  simultaneously 
with  the  optica]  imaging  measurements  seems  promising  for  validation  of  in  vivo  measurements.  In  summary,  time- 
resolved  spectroscopic  imaging  has  the  ability  to  detect  tumors  in  tissues  and  will  be  useful  in  3-D  image  reconstruction. 
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